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Introduction: 

-  Tumor  blood  volume  and  micro-vascular  density  are  parameters  anatomically  and  functionally 
associated  with  tumor  angiogenesis.  During  the  last  decade,  rigorous  modeling  of  the  light 
propagation  in  the  near  infrared  region,  combined  with  the  advancements  of  light  source  and 
detectors,  has  improved  the  diffused  light  measurements  and  made  possible  the  application  of 
tomographic  techniques  for  characterizing  and  imaging  tumor  angiogenesis  [1-6].  If  a  single 
wavelength  is  used,  optical  absorption  related  to  angiogenesis  and  other  normal  blood  vessels 
can  be  measured.  If  two  or  more  optical  wavelengths  are  used,  both  oxy-hemoglobin  and 
deoxy-hemoglobin  concentrations  can  be  measured  simultaneously.  In  addition,  optical 
scattering  is  related  to  cell  activations  because  scattering  is  sensitive  to  any  particles  that  are  of 
the  size  of  the  optical  wavelength.  However,  the  NIR  technique  has  not  been  widely  used  in 
clinics  and  the  fundamental  problem  remains  the  intense  light  scattering.  As  a  result,  diffused 
light  probes  a  wide  region  instead  of  propagating  along  a  straight  line.  Localization  or  imaging 
based  on  tomographic  inverse  scattering  approaches  suffers  from  low  spatial  resolution  and  the 
inversion  problem  is,  in  general,  underdetermined  and  ill-posed. 

Currently,  ultrasound  is  routinely  used  in  conjunction  with  mammography  to  differentiate  simple 
cysts  from  solid  lesions.  When  the  criteria  for  a  simple  cyst  are  strictly  adhered  to,  the  accuracy 
of  ultrasound  is  96%-100%  .  However,  the  ultrasound  appearance  of  benign  and  malignant 
lesions  has  considerable  overlapping  features  [7-8],  which  has  prompted  many  radiologists  to 
recommend  biopsies  on  most  solid  nodules.  This  results  in  a  large  number  of  biopsies  yielding 
normal  or  benign  breast  tissue  (currently  70%  to  80%  of  biopsies  are  normal  [9] ).  In  addition, 
the  diagnostic  accuracy  of  ultrasound  depends  largely  on  the  experience  of  physicians.  The  best 
results  reported  in  the  literature  on  distinguishing  benign  from  malignant  solid  breast  lesions 
were  given  by  Stavros  et  al  [10].  In  this  reference,  a  total  of  750  palpable  solid  breast  lesions 
were  studied.  Despite  the  known  overlap  features  in  some  lesions,  ultrasound  was  able  to 
correctly  classify  123  of  1 25  malignant  lesions  as  intermediate  or  malignant. 

We  have  introduced  the  use  of  optical  imaging  as  an  adjunct  to  ultrasound  in  differentiating 
benign  from  malignant  lesions  [1 1].  We  have  demonstrated  that  optical  contrast  between  benign 
and  malignant  lesions  can  provide  more  diagnostic  information  to  further  classify  the  acoustic 
intermediate  group,  thereby  improving  ultrasound  specificity  and  reducing  unnecessary  biopsies. 

We  have  introduced  a  novel  combined  imaging  using  lesion  structure  images  provided  by  co¬ 
registered  ultrasound  to  improve  NIR  functional  imaging  capability  in  determining  optical 
properties  [12].  With  the  a  priori  knowledge  of  lesion  location  and  shape,  NIR  imaging 
reconstruction  can  be  localized  within  specified  spatial  and  temporal  regions.  As  a  result,  the 
reconstruction  is  over-determined  because  the  total  number  of  unknown  optical  properties  is 
reduced  significantly.  In  addition,  the  reconstruction  is  less  sensitive  to  noise  because  the 
convergence  can  be  achieved  within  a  small  number  of  iterations.  We  have  conducted  a  series  of 
experiments  to  assess  the  improvement  on  reconstructed  optical  absorption  coefficients  with  the 
a  priori  target  depth  [13]  and  spatial  location  information  [14]. 

In  this  DOD  ARMY  sponsored  study,  we  proposed  to  achieve  the  flowing  objectives 

1)  to  refine  our  existing  NIR  optical  imaging  system  hardware  toward  high  signal-to-noise  ratio 

and  fast  data  acquisition  (task  1); 

2)  to  implement  to  imaging  software  for  clinical  studies  (task  1) 

3)  to  optimize  the  combined  probe  design  through  simulation  and  phantom  experiments  (task  1), 
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4)  to  validate  the  combined  imaging  in  cancer  detection  and  diagnosis  through  clinical  studies 
(task  2). 

We  have  successfully  completed  task  1  proposed  in  the  application  during  the  first  year  [13-14] 
and  have  partially  completed  task  2  by  recruiting  a  total  of  24  patients  and  correlating  the  NIR 
optical  imaging  results  with  ultrasound  findings  [15].  The  preliminary  results  are  summarized  in 
the  body  of  the  report.  More  patients  are  being  recruited  and  statistics  will  be  reported  shortly. 


Body 

Under  the  supports  of  DOD  ARMY  we  have  completed  the  combined  probe  and  prototype  NIR 
imager.  This  imager  consists  of  12  dual  wavelength  sources  and  8  optical  detectors  [14],  and 
the  sources  and  detectors  are  coupled  through  optical  fibers  to  a  hand-held  probe  shown  in 
Figure  1(a)  and  the  probe  dimensions  are  shown  in  Figure  1(b).  A  commercial  ultrasound  probe 
(7  MHz  linear  array)  is  simultaneously  deployed  in  the  middle  of  the  combined  probe. 
Currently,  this  probe  is  used  at  the  Cancer  Center  of  University  of  Connecticut  Health  Center  for 
the  proposed  clinical  studies.  With  this  combined  probe,  we  have  demonstrated  the  advantage  of 
dual-modality  diagnosis  [15]. 


Fig.l.  (a)  Picture  of  a  hand-held  combined  probe  and  a  frequency  domain  NIR  imager,  (b)  Sensor 
distribution  of  the  combined  probe.  The  diameter  of  the  combined  probe  is  10  cm.  Smaller  circles  are 
optical  source  fibers  and  big  circles  are  detector  fibers.  A  commercial  ultrasound  probe  is  located  at  die 
center  and  its  dimensions  are  5.6  cm  by  1  cm. 

Before  the  clinical  studies,  we  have  fully  tested  our  imaging  algorithms  with  phantom  targets 
which  emulate  tumor  acoustic  and  optical  properties.  A  new  method,  moment  method,  has  been 
invented  [16]  and  tested  with  phantom  targets.  However,  when  we  acquired  first  patient  data 
from  several  volunteers  in  November  2001  and  used  existing  algorithms  as  well.as  our  moment 
method  for  optical  imaging,  we  had  difficulties  in  obtaining  consistent  results  for  all  the  cases. 
We  further  investigated  the  imaging  algorithms  and  have  successfully  invented  a  two-step 
imaging  reconstruction  algorithm  reported  in  Ref.  [15].  This  algorithm  provides  consistent 
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•wavelength-dependent  optical  absorption  and  hemoglobin  concentration  of  a  large  breast 
'Carcinoma  with  necrotic  core.  To  the  best  of  our  knowledge,  such  detailed  distributions  have 
not  been  achieved  by  NIR  only  imaging  reported  in  the  literature. 


An  example  from  Ref.  15  is  given  here.  Figure  2  (a)  shows  a  gray  scale  ultrasound  image  of  a 
palpable  lump  of  a  44-year-old  woman.  The  lesion  was  located  at  6  to  8  o’clock  position  of  the 
left  breast  at  approximately  1.5  cm  depth.  Ultrasound  showed  an  irregular  poorly  defined 
hypoechOic  mass  and  the  lesion  was  considered  as  highly  suspicious  for  malignancy.  An 
ultrasound  guided  core  needle  biopsy  was  recommended.  Biopsy  results  revealed  that  the  lesion 
was  a  high  grade  in-situ  ductal  carcinoma  with  necrosis. 

Multiple  optical  measurements  at  two  orthogonal  positions  were  simultaneously  made  with 
ultrasound  images  at  the  lesion  location  as  well  as  at  approximately  the  same  location  of  the 
contralateral  normal  breast.  The  image  reconstruction  was  performed  using  the  NIR  data 
simultaneously  acquired  with  the  ultrasound  image  shown  in  Fig.2  (a).  The  optical  image 
reconstruction  revealed  optimal  lesion  centers  at  approximately  (-1.1  cm,  0.3  cm,  1.7  cm)  for  780 
nm  and  (-0.9  cm,  -0.7  cm,  1.7  cm)  for  830  nm,  respectively,  and  optimal  diameters  of  4.28  cm  x 
5.18  cm  x  1.96  cm.  The  detailed  absorption  maps  with  high  absorption  non-uniformly 
distributed  around  the  lesion  boundaries  at  both  wavelengths  are  shown  in  Fig.  2  (b)  and  (c).  By 
assuming  that  the  major  chromophores  are  oxygenated  (oxyHb)  and  deoxygenated  (deoxyHb) 
hemoglobin  molecules  in  the  wavelength  range  studied,  we  can  estimate  the  distribution  of  total 
hemoglobin  concentration  as  shown  in  Fig.2  (d).  The  extinction  coefficients  used  for 

calculating  oxyHb  and  deoxyHb  concentrations  were  e^0 =2.544,  £^°02=  1.695,  £Hb°=1.797, 

£Hb02=2.419  in  a  natural  logarithm  scale  with  units  of  inverse  millimoles  times  inverse 

centimeters.  The  measured  average  cancer  and  background  total  hemoglobin  concentrations 
were  116.9  p  moles  and  47.6  |i  moles,  respectively.  These  values  agree  with  3-D  average 
measurements  reported  in  the  literature. 


It  is  interesting  to  note  that  the  absorption  distributions  at  both  wavelengths  as  well  as  total 
hemoglobin  concentration  were  distributed  heterogeneously  at  the  cancer  periphery.  To  the  best 
of  our  knowledge,  such  fine  distributions  have  not  been  reported  by  using  NIR  only 
reconstruction  techniques.  However,  this  finding  agrees  with  the  published  literature  showing 
that  breast  cancers  have  higher  blood  volumes  than  non-malignant  tissue  due  to  angiogenesis, 
especially  at  the  cancer  periphery  [17].  In  addition,  the  carcinoma  reported  here  had  a  necrotic 
core  which  could  lead  to  low  absorptions  observed  at  both  wavelengths  in  the  center  region. 
Ideally,  pathologic  report  with  blood  vessel  density  count  could  be  used  to  correlate  the  NIR 
findings  reported  here.  Future  plans  are  to  obtain  microvessel  counts  in  the  region  of  the  tumor 
as  well  as  surrounding,  non-malignant  breast  tissue. 


Fig.  2.  (a)  A  gray  scale  ultrasound  image  of  a  palpable  lump  of  a  44-year-old  woman. 
Ultrasound  showed  an  irregular  poorly  defined  hypoechoic  mass  and  the  lesion  was  considered 
as  highly  suspicious  for  malignancy.  Reconstructed  optical  absorption  maps  at  780  nm  (b)  and 
830  nm  (c).  Vertical  color  bars  are  absorption  coefficient  in  die  unit  of  cm'1,  (d)  is  the  total 
hemoglobin  concentration  map.  The  vertical  color  bars  are  p.  moles.  The  NIR  data  were 
simultaneously  acquired  with  ultrasound  image  shown  in  Fig.2  (a).  Each  image  consists  of  7 
slices  obtained  in  0.5  cm  spacing  from  0.5  cm  to  3.5  cm  in  depth.  The  vertical  and  horizontal 
axes  correspond  to  x  and  y  dimensions  of  9  cm  by  9  cm. 
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In  summary,  we  have  reported  a  two-step  image  reconstruction  scheme  and  we  have 
demonstrated  that  aided  with  co-registered  ultrasound,  detailed  optical  absorption  and  total 
hemoglobin  distributions  of  the  breast  carcinoma  can  be  obtained.  These  distributions  reveal 
tumor  angeogenesis  and  provide  valuable  information  for  breast  cancer  diagnosis  and  treatment. 

We  have  studied-  a  total  of  24  patients  for  the  past  year  and  are  still  recruiting  more  patients  to 
our  study.  Patient  recruiting  is  more  difficulty  than  we  have  anticipated.  First,  we  have  to  study 
patients  who  have  suspicious  lesions  seen  by  mammograms  and/or  ultrasound.  Second,  some  of 
the  good  candidates  may  not  want  to  participate,  particularly  cancer  patients.  We  have  received 
one-year  extension  on  this  contract  and  we  hope  to  recruit  as  many  patients  as  possible  in  the 
coming  year  and  to  establish  statistic  significance  of  our  combined  imaging  method. 

Reportable  outcomes 

Two  important  issues  directly  related  to  combined  breast  imaging  are  probe  geometry  and  optical 
data  acquisition  system. 

Both  reflection  and  transmission  geometries  have  been  used  by  researchers  for  optical 
tomography  and  .transmission  geometry  is  more  popular  than  reflection  due  to  dynamic  range 
considerations.  In  reflection  geometiy,  patients  are  placed  in  supine  position  and  the  imaging 
probe  consisting  of  multiple  source  fibers  and  detector  fibers  is  placed  on  top  of  the  breast. 
Reflected  light  detected  from  multiple  detector  locations  is  used  for  image  reconstruction. 
However,  the  dynamic  range  of  detected  signals  is  very  large  because  lesions  could  be  closer  to 
the  skin  or  deeply  seated.  We  have  chosen  reflection  geometry  in  our  current  hand-held  probe 
configuration  because  conventional  ultrasound  probes  are  used  in  reflection  mode.  Therefore, 
our  combined  probe  design  compromises  both  ultrasound  tomography  and  optical  tomography, 
hi  transmission  geometry,  patients  are  in  prone  position  and  the  examined  breast  is  inserted  in  a 
liquid  reservoir.  Optical  source  and  detector  fibers  can  be  distributed  around  the  liquid  reservoir. 
Due  to  gravity,  the  breast  is  much  longer  in  the  direction  perpendicular  to  the  source-detector 
plane.  As  a  consequence,  the  depth  of  the  lesion  relative  to  the  source  and  detector  fibers  does 
not  change  dramatically  and  the  dynamic  range  of  the  detected  signals  is  reduced.  However, 
ultrasound  probes  are  not  commonly  used  in  transmission  geometry  due  to  extra  distortion  of 
acoustic  wave  between  the  liquid  and  the  breast.  We  have  devoted  part  of  our  efforts  to  optical 
•  transmission  geometry  study  in  an  attempt  to  find  out  the  best  compromise  design. 

Our  results  are  obtained  with  2-dimensional  transmission  geometry  set-up  and  have  showed  that 
similar  performance  can  be  expected  as  long  as  a  priori  information  of  lesion  locations  is  used  in 
optical  image  reconstruction  [18]. 

Regarding  optical  data  acquisition  system,  frequency  domain  and  DC  systems  are  frequently  used 
by  many  research  groups.  The  advantage  of  using  frequency  domain  system  is  that  both 
amplitude  and  phase  information  can  be  obtained  and,  therefore,  simultaneous  reconstruction  of 
absorption  and  scattering  coefficients  is  feasible.  The  disadvantage  is  the  complexity  in 
constructing  high  frequency  transmission  and  reception  circuitries,  particularly  when  many 
parallel  channels  are  needed  for  fast  data  acquisition.  In  addition,  the  system  cost  is  high.  DC 
systems  have  advantages  of  simplicity  in  circuit  construction  and  low  cost.  The  drawbacks  are 
that  the  phase  information  is  not  available  and  the  cross-talk  between  the  reconstructed  absorption 
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'and  scattering  coefficients  may  be  unavoidable.  However,  when  reconstruction  of  lesion 
'absorption  maps  is  concerned,  frequency  domain  and  DC  systems  may  demonstrate  similar 
performance.  We  have  compared  the  performance  of  both  systems,  in  constructing  absorption 
maps  under  the  conditions  of  NIR  reconstruction  only  and  ultrasound  guided  NIR  reconstruction 
and  the  results  as  we  expected  are  reported  in  [18]. 


Key  Research  Accomplishments: 

•  invented  a  new  optical  imaging  algorithm  which  can  be  used  successfully  for  in  vivo  data 

•  acquired  data  from  24  patients  with  suspicious  lesions 

•  studied  combined  imaging  using  transmission  geometry  in  an  attempt  to  find  out  the  best 
compromise  in  the  imaging  probe  design. 


Conclusions: 

We  have  simultaneously  acquired  ultrasound  and  NIR  optical  data  from  24  patients  and  have 
showed  encouraging  results  for  imaging  breast  lesion  agiogenesis  using  the  combined  approach. 
In  addition,  we  have  studied  transmission  geometry  for  combined  imaging  in  an  attempt  to  find 
out  the  best  geometry  for  using  the  combined  approach. 
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NIR  imaging  reconstruction  with  ultrasound  localization:  feasibility  study  using  transmission 

geometry 

Minming  Huang,  Tuqiang  Xie,  Nan  Guang  Chen  and  Quing  Zhu* 

Electrical  and  Computer  Engineering  Department 
University  of  Connecticut,  Storrs,  CT  06269 

Abstract 

In  this  paper  we  report  experimental  results  of  NIR  imaging  reconstruction  with  ultrasound 
localization.  NIR  data  were  obtained  from  frequency  domain  and  DC  systems  with  source  and 
detector  fibers  configured  in  transmission  geometry.  Both  high  and  low  contrast  targets  located 
closer  to  the  boundary  and  closer  to  the  center  of  the  turbid  medium  were  reconstructed  using 
NIR  rdata  only  and  NIR  data  with  ultrasound  localization.  Results  have  shown  that  the 
improvements  of  reconstructed  mean  absorption  coefficients  of  high-contrast  targets  range  from 
20-26%  to  67-87%  with  ultrasound  localization.  The  measured  FWHMs  have  been  improved 
from  210-250%  to  110-120%  of  the  true  target  size  with  ultrasound  localization.  For  low 
contrast  target  cases,  the  targets  were  hardly  recognizable  for  both  on  center  and  off-center  cases 
with  NIR  data  obtained  from  both  systems.  With  ultrasound  localization,  the  reconstructed  mean 
absorption  coefficients  were  within  67%  of  the  true  value  and  the  target  can  be  visualized.  When 
the  change  of  target  scattering  coefficient  was  small  from  the  background  medium,  the 
performance  of  the  DC  system  was  comparable  to  that  of  the  frequency  domain  system. 
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Introduction  9 

Functional  imaging  with  near  infrared  (NIR)  light  has  found  potential  applications  in  many  areas. 
u9  However,  one  fundamental  problem  that  we  have  to  overcome  is  the  intense  scattering  of 
light.  As  a  result,  diffusive  light  probes  a  widespread  region  instead  of  providing  information 
along  a  straight  line.  Multiple  measurements  are  always  correlated  as  a  result  of  the  overlapping 
of  probed  regions.  Therefore,  increasing  the  total  number  of  measurements  does  not  necessarily 
provide  more  independent  information  for  image  reconstruction.  In  general,  the  inverse  image 
reconstruction  is  underdetermined  and  ill-posed.  The  behavior  of  reconstruction  algorithms  is 
affected  by  many  factors,  such  as  system  signal  to  noise  ratio,  probe  configuration,  and  the 
regulation  schemes  used  in  image  reconstructions. 

Current  approaches  of  image  reconstruction  algorithms  are  1)  simple  back  projection  methods,34 
2)  perturbation  methods,10'12  and  3)  finite  element  methods  (FEM).13'15  The  back  projection 
method  provides  real  time  estimation  of  coarse  optical  properties  of  lesions.  However,  the 
reconstructed  image  resolution  was  low  and  the  lesions  that  appeared  in  images  were  often 
displaced  from  true  locations.  Perturbation  methods  are,  in  general,  based  on  linear 
approximations  to  the  heterogeneous  functions.  Bom  and  Rytov  approximations  are  examples. 
Measurements  made  between  the  background  and  the  heterogeneous  medium  were  used  to  relate 
optical  signals  ait  the  measurement  surface  to  absorption  and  scattering  variations  in  each  volume 
element  within  the  sample.  The  least  square  method  was,  in  general,  used  to  formulate  the 
inverse  problem.  Iterative  search  methods,  such  as  conjugate  gradient  techniques,  were 
employed  to  iteratively  solve  the  inverse  problem.  However,  accurate  estimation  of  the  target 
optical  properties  depends  on  the  accurate  estimation  of  the  background  optical  properties,  which 
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are,  in  general,  not  easy  to  obtain  in  breast  tissues.  In  addition,  when  the  absorption  or  the 
scattering  coefficients  of  the  lesions  are  significantly  higher  than  the  background,  the  linear 
perturbation  methods  will  not  give  accurate  optical  properties.  Reconstructions  based  on  finite 
element  methods  provide  higher  order  estimations  to  heterogeneous  functions.  However,  as  the 
scattering  and  absorption  coefficients  are  expanded  over  local  basis  functions,  the  number  of 
unknowns  is  increased  considerably.  Nonetheless,  the  independent  information  from  diffuse 
light  is  limited,  and  it  will  essentially  stop  to  increase  when  the  source-detector  pairs  reach  a 
certain  number.  In  this  case,  the  inverse  problem  will  become  more  ill-posed,  and  the  behavior 
of  the  reconstruction  algorithm  may  be  unpredictable. 

Reconstructions  with  the  aid  of  a  priori  target  geometry  information  provided  by  co-registered 
ultrasound  have  shown  promising  results  in  improving  the  accuracy  of  reconstructed  optical 
properties  and  the  localization  of  targets.1619  In  this  method,  an  ultrasound  probe  and  NIR 
source  and  detector  fibers  have  been  deployed  on  a  hand-held  probe  and  configured  in  the 
reflection  geometry.  The  a  priori  tissue  type  and  lesion  location  as  well  as  shape  provided  by 
co-registered  ultrasound  can  guide  the  NIR  image  reconstruction  to  localized  target  regions. 
Therefore,  the  total  number  of  voxels  with  unknown  optical  properties  can  be  reduced 
significantly  and  the  inversion  becomes  overdetermined.  In  general,  the  solution  is  unique  and 
the  iterative  algorithms  converge  very  fast. 

In  this  paper  we  explore  experimentally  the  utility  of  the  combined  imaging  with  NIR  sources 
and  detectors  configured  in  transmission  geometry.  The  standard  pulse-echo  technique  is  used 
to  obtain  co-registered  ultrasound  images.  Compared  with  NIR  reflection  geometry. 
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transmission  geometry  may  have  the  advantage  of  reduced  target  dynamic  range  and  higher 
sensitivity.  In  the  reported  experiments,  the  FEM  forward  solver  is  used  to  generate  a  Jacobian 
weight  matrix  for  the  NIR  image  reconstruction,  and  the  inversion  is  performed  by  using  NIR 
data  only  and  NIR  data  with  a  priori  target  geometry  provided  by  co-registered  ultrasound. 

Frequency  domain  and  DC  systems  are  frequently  used  by  many  research  groups.  The 
advantage  of  using  frequency  domain  system  is  that  both  amplitude  and  phase  information  can 
be  obtained  and,  therefore,  simultaneous  reconstruction  of  absorption  and  scattering  coefficients 
is  feasible.  The  disadvantage  is  the  complexity  in  constructing  high  frequency  transmission  and 
reception  circuitries,  particularly  when  many  parallel  channels  are  needed  for  fast  data 
acquisition.  This  situation  is  more  pronounced  when  nearly  100  source  and  detector  positions 
are  needed  for  3-D  imaging  in  transmission  geometry.  In  addition,  the  system  cost  is  high.  DC 
systems  have  advantages  of  simplicity  in  circuit  construction  and  low  cost.  The  drawbacks  are 
that  the  phase  information  is  not  available  and  the  cross-talk  between  the  reconstructed 
absorption  and  scattering  coefficients  may  be  unavoidable.  However,  where  reconstruction  of 
lesion  absorption  maps  is  concerned,  frequency  domain  and  DC  systems  may  demonstrate 
similar  performance.  In  this  paper,  we  have  compared  the  performance  of  both  systems  in 
constructing  absorption  maps  under  the  conditions  of  NIR  reconstruction  only  and  ultrasound 
guided  NIR  reconstruction. 


2.  Forward  model 


We  have  adopted  a  commercial  CFD  package  FLUENT,  which  allows  users  to  change  the 
governing  equation  through  user  defined  source  functions.  The  steady-state  governing  equation 
for  transport  of  a  scalar  quantity  <p  in  FLUENT  is 

-^-(rR|^-)  =  S  k  =  l . N  (1) 

OX;  OX;  * 

where  TR  is  the  diffusion  coefficient  of  <j> ,  Sf  is  the  source  of  <0  per  unit  volume,  and  k  is  the 
index  of  image  voxel  number.  We  have  replaced  the  source  in  Eq.  (1)  by  the  user  defined 

A 

source  S&  =  S#  -  jl$k  and  obtained  a  new  governing  equation  as 

d  <%■  A 

-^7(D^)+^k=s-  (2) 

1  1 

This  is  the  photon  diffusion  equation,  where  //a  is  the  absorption  coefficient  and  TR  is  replaced 
by  the  photon  diffusion  coefficient  D  =  — - - .  it’  is  the  reduced  scattering  coefficient. 

3(Ht+M 


To  obtain  the  frequency  domain  forward  model,  we  have  used  the  transient  governing  equation 


at 


0X; 


)+pA=s^ 


to  obtain  time  domain  waveforms  and  then  calculate  the  DC  component,  amplitude  and  phase  of 
each  waveform.  The  step  used  in  the  time  domain  waveform  calculation  is  1/3  of  the  period  T, 
where  T  =  1/f  and  f  =  140MHz,  which  is  the  modulation  frequency  of  our  experimental  system. 
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Assume  that  the  change  of  D  is  small  and  the  scattering  wavers  primarily  generated  by  the 

A<J> 

absorption  inhomogeneity.  The  Jacobian  matrix  Wy  =  — —  that  relates  the  scattering  wave  at 

Apaj  ’ 

boundary  cell  i  and  imaging  voxel  j  with  absorption  coefficient  change  A|i.aj ,  is  given  as 


3.  Computation  procedures 

We  have  used  a  two-dimensional  mesh  with  4106  triangle  elements  and  2130  nodes  (see  Fig  1). 
The  radius  of  physical  boundary  is  chosen  as  45.72  cm  and  the  extrapolation  distance 

(L  =  — ^ — )  is  calculated  as  0.167  cm,  where  U,+Ua  is  the  background  value  and  is  chosen 
as  6  cm 


In  Jacobian  matrix  calculation,  we  have  evaluated  the  linearity  of 


A^aj 


verse  A(Xaj  with 


paj  changing  from  0.02  cm'1  to  0.22  cm'1  at  many  cell  locations.  The  linearity  is  related  to  the 
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distance  between  the  cell  and  the  source,  and  the  linearity  is  better  when  the  distance  is  larger. 
The  worst  case  is  that  the  cell  j  is  next  to  the  source  cell.  For  the  DC  component,  this 
dependence  is  not  significant,  and  the  error  due  to  the  linearity  approximation  is  always  less  than 
0.5%.  For  frequency  domain  calculations,  the  real  and  imaginary  part  of  of  several  cells  are 

one  or  two  orders  smaller,  and  they  are  very  sensitive  to  small  perturbations.  However,  even  for 
the  worst  case,  the  error  due  to  the  use  of  linear  approximation  is  also  less  than  0.7%  for  306  out 
of  308  cells  at  source-detector  layer  (See  Fig  1).  Therefore,  we  have  simplified  the  Jacobian 

matrix  calculation  by  computing  — —  at  two  points  0.02  cm'1  and  (J.^  =0.22  cm '.  Figure 

Apaj 

2  shows  verse  p.2J  with  paj  varies  from  0.02cm'1  to  paj=  0.22cm'1,  where  cell  j  is  located  at  the 
center.  Figure  2(a)  is  the  real  part  of  $g  verse  p.aj ,  and  Fig.  2  (b)  is  the  imaginary  part  of 
verse  p.aj .  Figure  2(c)  is  the  DC  component  of  4>tJ  verse  paj . 

The  experimental  data  have  included  many  unknown  systematic  factors  (unknown  source  and 
detector  gains,  etc.).  To  remove  those  unknown  factors,  we  have  used  the  normalized  difference 
method  (NDM)  given  in  Ref  20  as 

=  K!  Ew<.«)]  <5> 

where  4m0)  is  the  measured  heterogeneous  data  associated  with  the  source-detector  pair  i  with 
the  target  in  the  homogeneous  medium,  <j)mr(i)  is  the  measured  homogeneous  data,  <pcrQ)  is  the 
calculated  forward  homogeneous  data  used  for  weight  matrix  calculation.  0  is  a  complex  value 
in  frequency  domain  case.  The  unknown  system  factors  present  in  both  sets  of  measured 


*«(n  ~ 


mr  (i) 


cr  ( I ) 


mr  (i) 
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heterogeneous  and  homogeneous  data  are  cancelled  by  taking, the  ratio  of  the  perturbation 
(heterogeneous  -  homogeneous)  and  the  reference  (homogeneous)  measurements. 


The  conjugate  gradient  method  is  used  for  the  inversion.  Since  the  error  function 


n  w  4p 


V)  fmrjj)  2  reduces  consistently  and  becomes  flat  after  certain  iterations  (around 


a  0  cKO 

Ymrii) 


30  for  reconstruction  with  NIR  data  only,  and  3  for  NIR  with  ultrasound  guidance),  the  stopping 
criteria  are  chosen  to  be  100  and  3  iterations  for  NIR  only  and  NIR  with  ultrasound  localization, 
respectively.  The  corresponding  computation  time  is  around  3  minutes  for  NIR  data  only  and  30 
seconds  for  NIR  with  ultrasound  information. 


In  experiments,  data  are  obtained  from  the  3D  model  (finite  cylindrical  medium  and  finite  length 
target).  However,  the  Jacobian  matrix  computation  is  based  on  the  2D  model  (infinite  cylindrical 
medium  and  infinite  length  target).  According  to  Refs  21  and  22,  the  2D/3D  difference  between 
intensity  is  reasonably  independent  of  source-detection  separation,  and  it  is  reasonable  to  correct 
it  by  simply  multiplying  or  dividing  a  constant.  Therefore,  NDM  not  only  cancels  the  unknown 
source  detector  gains  but  also  partially  eliminates  the  2D/3D  model  mismatch  by  taking  the  ratio 
of  the  perturbation  and  the  reference  measurements. 

According  to  Refs  21  and  22,  the  2D/3D  model  mismatch  can  cause  some  artifacts  near  the 
boundary.  We  have  used  constraints  based  on  the  assumption  that  there  are  no  perturbations 
outside  the  detector  and  source  layer  to  partially  eliminate  the  artifacts. 
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4.  Experimental  systems  > 

To  compare  the  performance  of  frequency  domain  and  DC  systems  without  and  with  ultrasound 
localization,  we  have  constructed  both  systems.  The  details  of  our  multi-source  and  multi¬ 
detector  frequency  domain  system  can  be  found  in  Ref.  18.  Briefly,  it  consists  of  12  pairs  of 
dual  wavelength  sources  (780  nm  and  830  nm)  and  8  parallel  detectors.  The  light  sources  are 
laser  diodes  amplitude  modulated  at  140  MHz,  and  the  detectors  are  photon  multiplier  tubes 
(PMT’s).  These  sources  and  detectors  are  coupled  to  the  medium  through  optical  fibers. 

Our  new  DC  system  consists  of  6  pairs  of  dual-wavelength  laser  diodes  (780  nm  and  830  nm) 
amplitude  modulated  at  20  KHz  to  avoid  DC  fluctuations  (see  Fig.  3).  12  detector  fibers  are 
sequentially  coupled  to  a  PMT  through  a  mechanical  multiplexer,  which  uses  a  stepper  motor  to 
accurately  collimate  a  coupling  light  guide  on  a  rotation  plate  to  one  of  detection  fibers.  Each 
light  source  has  two  output  levels  of  30  dB  difference  to  avoid  PMT  saturation  and  improve  the 
system  dynamic  range.  A  National  Instrumentation  data  acquisition  (DAQ)  card  is  used  to 
generate  sinusoidal  waveforms  on  top  of  a  constant  current  to  drive  one  laser  diode  at  a  time. 
The  photon  density  wave  detected  by  the  PMT  is  amplified  and  then  digitized  by  the  same  DAQ 
card.  The  corresponding  amplitude  is  retrieved  from  the  acquired  waveform  on  a  PC.  The  digital 
ports  of  the  DAQ  card  are  used  to  switch  laser  diodes  sequentially,  control  the  output  level,  and 
input  a  feedback  signal  for  close-loop  stepper  motor  control.  The  total  data  acquisition  time  for  a 
complete  data  set  is  around  1  minute. 

Since  our  DC  system  is  modulated  at  20  kHz,  the  measured  amplitude  is  an  AC  signal  instead  of 
the  DC  component  we  used  for  forward  Jacobian  matrix  calculation.  We  have  estimated  the 
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difference  using  a  spherical  wave  of  light  energy  density  (e°kp*/Dp)  of  two  wave  numbers 
k,  =A/-pa  /D  andk2  =  ^/(— jxa  +  jx©/c)/D ,  where  co -In  -2xl04Hz,  c  and  D  are  the  speed 
of  light  and  the  reduced  diffusion  coefficient,  respectively.  The  normalized  difference 


p)  _  e(fap) 


can  be  approximated  as  ,  and  it  is  larger  when  the  wave 


„</M) 


2c 


propagates  deeper  into  the  tissue  and  the  background  absorption  is  smaller.  For  the  worst  case  of 
fjLa  =  0.01  cm'1  and  p  =  10cm,  the  difference  is  less  than  3.295E-5.  Therefore,  we  directly  use 

the  Jacobian  matrix  calculated  from  the  DC  component  to  relate  the  measured  amplitude  data  for 
the  inversion. 


During  experiments,  source  and  detector  fibers  were  placed  on  a  circular  plane  (see  Fig  4).  A 
commercial  ultrasound  probe  of  3.5  MHz  central  frequency  was  placed  on  the  top  of  a  water  tank 
filled  with  Intralipid.  This  one-dimensional  commercial  probe  provides  cross-section  image  in 
y-z  imaging  plane  (called  B-scan),  where  y  is  the  lateral  direction  and  z  is  the  propagation 
direction.  By  translating  the  probe  mechanically  in  x  direction,  we  acquired  the  3-D  volumetric 
image  data.  Windowing  the  3D  data  in  z  direction  at  a  particular  depth  provides  the  2D  target 
spatial  images  at  x-y  plane  (called  C-scan).  C-scans  are  co-registered  with  NIR  images.  0.6% 
Intralipid  solution  was  used  as  a  homogeneous  background.  The  fitted  na  and  //’  of  the 
Intralipid  at  780  nm  were  0.02cm'1  and  6cm'1,  respectively. 

Since  we  have  a  limited  number  of  source  and  detector  positions,  we  have  placed  8  source  fibers 
at  a  quarter  circle  and  8  detector  fibers  at  the  opposite  quarter  for  frequency  domain 
measurements  (see  Fig.5  (a)).  For  DC  measurements,  we  placed  6  source  fibers  at  a  quarter  circle 
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and  12  detector  fibers  at  the  opposite  quarter  (see  Fig.5  (b)).  After  obtaining  one  set  of 
heterogeneous  NIR  data  with  a  target  in  the  Intralipid  (8x8  measurements  for  the  frequency 
system,  6  x  12  measurements  for  the  DC  system),  we  simply  rotated  the  target  by  90°  for 
frequency  domain  system  and  77.140  for  the  DC  system  to  get  another  set  of  data,  respectively. 
Figure5  (c)  shows  the  rotation  scheme  used  during  frequency  domain  experiments.  Three 
rotations  were  performed  to  obtain  one  complete  data  set  for  reconstruction,  which  covered  360°. 
For  DC  measurements,  a  similar  rotation  scheme  was  used  and  four  rotations  were  performed  to 
obtain  one  complete  data  set.  3-D  ultrasound  data  were  obtained  simultaneously  with  NIR  data 
acquisition. 

5.  Results 

For  NIR  reconstruction  using  transmission  geometry,  the  target  dynamic  range  is  related  to  the 
location  of  the  target  and  the  contrast  of  the  target.  It  is  much  harder  to  reconstruct  a  low 
contrast  target  located  at  the  middle  of  the  2D  mesh.  In  our  experiments,  a  nearly  circular  target 
(high/  low  contrast)  with  1.0cm  diameter  was  placed  at  two  typical  locations  (at  the  center,  closer 
to  the  boundary).  The  target  was  made  of  an  acrylamide  gel 16  with  calibrated  na .  For  the  high 
contrast  case,  the  absorption  coefficient  of  the  target  was  0.20  cm'1,  which  was  0.18cm'1  above 
the  background.  For  the  low  contrast  case,  the  calibrated  absorption  coefficient  of  the  target  was 
0.10  cm'1,  which  was  0.08cm'1  above  the  background. 

A.  Reconstruction  using  frequency  domain  system 

Figure  6  shows  experimental  results  of  a  1.0  cm-diameter  cylindrical  target  of  fj.a=  0.2  cm'1 
located  closer  to  the  boundary.  The  reduced  scattering  coefficient  fls  ’  was  controlled  the  in  the 
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same  way  as  the  background  which  was  6  cm’1.  Figure  6  (a)  is  the  p.a  distribution  reconstructed 
using  NIR  data  only.  We  calculated  the  mean  value  for  the  pixels  with  absorption  coefficients 
greater  than  the  half  maximum  in  the  target  mass  region.  The  reconstructed  mean  absorption 
coefficient  was  only  24%  of  the  true  value,  and  the  full  width  at  half  maximum  (FWHM)  was 
around  2.2cm,  which  was  220%  of  the  true  target  size.  Figure  6  (c)  is  a  B-scan  ultrasound 
image,  and  Figure  6  (d)  is  the  C-scan  ultrasound  image  obtained  at  the  source-detector  plane. 
The  ultrasound  image  was  segmented  first  and  then  a  threshold  was  used  to  map  out  target  and 
non-target  regions.  The  inverse  reconstruction  was  localized  to  the  target  region.  Figure  6  (b)  is 
the  reconstructed  fia  distribution  with  co-registered  ultrasound  localization.  With  ultrasound 

localization,  the  accuracy  of  reconstructed  mean  fitt  was  improved  from  24%  to  75%  of  the  true 
vale.  The  FWHM  was  improved  to  1.16  cm,  which  was  only  16%  larger  than  the  true  value. 
Figure  6  (e)  shows  the  comparison  of  the  reconstructed  fia  distribution  along  the  x-axis  as 

indicated  by  vertical  arrows  in  (a),  and  Fig.  6  (f)  shows  the  comparison  of  reconstructed  fia 
distribution  along  the  y-axis  as  indicated  by  horizontal  arrows  in  (a). 

Figure  7(a)  shows  the  reconstructed  Ha  distribution  of  the  same  high  contrast  target  located  at 
the  middle.  As  one  can  see,  the  spatial  distribution  is  poor  and  the  reconstructed  mean  [la  was 
about  26%  of  the  true-  value.  Figure  7(b)  is  the  co-registered  C-scan  image  of  ultrasound  and 
Fig.7  (c)  is  the  reconstructed  NIR  absorption  map  with  ultrasound  localization.  The 
reconstructed  mean  fia  value  has  improved  to  26%  and  67%  of  the  true  value. 
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For  a  lower  contrast  target  ( fia  =  0. 1  cm'1)  located  at  different  locations,  the  reconstructed  images 

were  dominated  by  artifacts.  Figure  8  shows  the  reconstructed  absorption  maps  without  (part 
(a))  and  with  (part  (b))  ultrasound  localization  when  the  target  was  located  closer  to  the 
boundary.  Without  ultrasound,  the  target  cannot  be  recognized.  With  ultrasound,  the  mean 
reconstructed  absoiption  coefficient  was  78%  of  the  true  value. 

Table  I  summaries  the  experimental  results  of  two  cases  when  the  target  was  located  at  different 
positions.  In  general,  when  the  target  was  located  at  the  middle,  the  results  obtained  with  NIR 
data  only  were  slightly  poorer  than  those  when  the  target  was  located  closer  to  the  boundary,  and 
the  improvements  with  ultrasound  localization  were  similar  to  those  when  the  target  was  located 
at  the  boundary. 

B.  Reconstruction  using  the  DC  system 

Figure  9  shows  experimental  results  of  the  1.0  cm-diameter  cylindrical  target  of  jua=  0.2  cm'1 
located  closer  to  the  boundary.  Figure  9  (a)  is  the  na  distribution  reconstructed  using  NIR  data 

only.  The  reconstructed  mean  absorption  coefficient  was  22%  of  the  true  value,  and  the  FWHM 
was  around  2.5cm,  which  was  2.5  times  that  of  the  true  target  size.  Figure  9  (b)  is  the  co¬ 
registered  C-scan  ultrasound  image  obtained  at  the  source-detector  plane.  Figure  9  (c)  is  the 
reconstructed  na  distribution  with  ultrasound  localization.  The  accuracy  of  reconstructed  mean 

Ha  was  improved  to  77%  of  the  true  value.  The  measured  FWHM  with  ultrasound  localization 
was  1.18  cm,  which  was  only  18%  larger  than  the  true  value.  Figure  9  (d)  and  (e)  show  the 
comparison  of  the  reconstructed  fia  distributions  along  x-axis  and  y-axis,  respectively. 
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The  reconstructed  absorption  map  is  dominated  by  artifacts  for*the  lower  fia  (0. 1  cm'1)  case. 
Figure  10  shows  a  low  contrast  target  located  closer  to  the  boundary.  The  target  can  only  be 
visualized  with  ultrasound  localization. 

Table  II  summaries  experimental  results  of  high  and  low  contrast  cases  when  the  target  was 
located  at  different  positions.  The  reconstructed  absorption  coefficients  and  improvements  with 
ultrasound  localization  are  comparable  to  those  obtained  from  the  frequency  domain  system. 

Comparing  results  obtained  from  frequency  domain  and  DC  systems,  we  find  no  significant 
difference  with  respect  to  accuracy  and  target  spatial  distributions.  The  reconstructed  absorption 
values  of  using  NIR  data  only  obtained  from  both  systems  are  very  close.  With  ultrasound 
localization,  the  mean  absorption  coefficients  obtained  by  the  DC  system  are  about  9-13  % 
higher  than  the  frequency  domain  results  for  the  low  contrast  target  case.  For  the  high  contrast 
target  case,  the  reconstruction  results  are  similar. 

6.  Discussions 

Based  on  our  experimental  results,  the  reconstructed  values  by  using  NIR  only  are  certainly  low 
and  the  target  spatial  distributions  are  poor.  These  problems  can  be  partially  improved  by 
iteratively  updating  the  Jacobian  matrix  to  account  for  higher  order  terms  of  the  weight  matrix. 
However,  this  procedure  is  time-consuming  and  prohibits  near  real-time  NIR  image  processing. 
With  the  guidance  of  spatial  target  distribution  provided  by  ultrasound,  the  reported  mean 
reconstruction  results,  are  within  67%  to  87%  .  of  the  true  absorption  coefficients  when  the  first 
order  linear  approximation  is  used.  Iterative  updating  of  the  localized  Jacobian  matrix  could 
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further  improve  the  accuracy  of  reconstructed  absorption  coefficients  with  reduced  computation 
load.  Currently,  we  are  pursuing  this  research  and  results  will  be  reported  in  the  near  future. 


In  the  Jacobian  matrix  computation,  we  have  calculated  — —  cell  by  cell  by  assuming  that  high 

^M’aj 

order  cross-talk  terms  between  cells  are  negligible.  This  assumption  is  equivalent  to  using  a 
linear  approximation  of  A0  as 
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where  index  i,  j  . 1  are  cell  index  numbers.  To  estimate  the  error  due  to  the  linear 

approximation,  we  calculated  the  difference  between  A<j)  and  its  linear  approximation  at  several 
locations  (next  to  the  source,  1cm  away  from  the  source,  at  the  middle  of  the  mesh).  The  error  is 
less  than  1 .8%  for  four  cells  (0.3  cm  target  size),  7.9%  for  eight  cells  (0.6  cm  target  size),  and 
27%  for  28  cells  (1 .0  cm  target  size). 


In  the  reported  phantom  studies,  we  assume  that  the  lesions  are  isolated  and  are  imbedded  in  a 
homogeneous  background.  Therefore,  we  localize  the  image  reconstruction  by  using  target 
geometry  obtained  from  ultrasound.  In  the  clinical  breast  imaging  cases,  the  background  tissues 
also  scatter  and  absorb  the  diffusive  light,  and  constraints  on  the  background  perturbation  will  be 
used  to  account  for  the  scattering  and  absorption  contributions  from  the  background.* 


In  this  paper,  the  combined  imaging  has  been  demonstrated  using  transmission  geometry  for 
imaging  absorbers.  In  principle,  the  combined  imaging  is  possible  for  imaging  scattering  targets 
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as  well.  Simultaneous  reconstruction  of  absorbing  and  scattering  contrasts  are  more  difficult 
because  of  the  crosstalk  between  them,  particularly  when  the  DC  system  is  used.  However, 
further  investigation  will  be  conducted  along  this  direction. 

Due  to  the  fact  that  the  2D/3D  data  difference  is  not  exactly  linear  dependent,  some  artifacts  are 
visible  in  the  boundary  regions  of  reconstructed  absorption  maps  when  NIR  data  are  used.  For 
the  lower  contrast  target  case,  this  effect  is  dominant  and  the  target  can  only  be  reconstructed  by 
using  a  priori  target  information  provided  by  ultrasound.  Since  photons  travel  in  three 
dimensions,  3D  meshes  are  more  appropriate  and  we  are  currently  working  on  3-D  model 
constructions  and  experimentations. 

7.  Summary 

This  paper  reports  experimental  results  of  NIR  imaging  reconstruction  with  ultrasound 
localization.  Transmission  geometry  was  used  for  NIR  imaging  reconstruction  and-standard 
pulse-echo  ultrasound  was  used  for  obtaining  co-registered  ultrasound  images.  The  results  have 
shown  that  reconstructed  mean  absorption  coefficients  of  phantom  targets  have  improved  from 
24%  to  75%  and  26%  to  67%  for  a  high-contrast  target  located  closer  to  the  boundary  and  closer 
to  the  center  of  the  turbid  medium,'  respectively,  when  a  frequency  domain  system  is  used. 
Similarly,  the  reconstructed  mean  absorption  coefficients  have  improved  from  22%  to  77%  and 
20%  to  76%  for  the  same  high-contrast  target  located  closer  to  the  boundary  and  closer  to  the 
center  of  the  turbid  medium,  respectively,  when  a  DC  system  is  used.  The  reconstructed 
absorption  maps  were  poor  in  general  when  targets  were  located  at  the  center  of  the  turbid 
medium.  Ultrasound  localization  improves  the  accuracy  of  reconstructed  values  and  reduces 
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image  artifacts.  With  ultrasound  localization,  the  measured  FWHMs  have  been  improved  from 
210-250%  to  1 10-120%  of  the  true  target  size.  For  a  low  contrast  target,  the  target  was  hardly 
recognizable  for  both  closer  to  boundary  and  closer  to  center  cases  with  NIR  data  obtained  from 
both  systems.  With  ultrasound  localization,  the  reconstructed  absorption  coefficients  were 
within  67%  of  the  true  value  and  the  target  can  be  visualized.  With  ultrasound  localization,  the 
reconstruction  speed  has  improved  by  a  factor  of  10  and  near  real-time  optical  imaging  becomes 
feasible.  Under  the  condition  that  the  change  of  scattering  coefficient  is  small,  the  performance 
of  the  DC  system  is  comparable  to  that  of  the  frequency  domain  system. 
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Figure  Captions 

Figure  1 .  Two  dimensional  mesh. 

Figure  2  versus  Ap,  with  p,,  changing  from  0.92  cm  '  .0  0.22  cm  '  for  a  cel,  locared  in  the 
AM-,j 

middle  of  die  mesh,  (a)  Real  par.  of  *  versus  p„ .  0»  Imaginary  par.  of  versus  P, . 

(c)  Calculation  is  made  using  amplitude  of  DC  data. 


Figure  3.  Schematic  of  our  new  DC  system. 

AO:  analog  output;  AI:  analog  input;  DO:  digital  output;  DI:  digital  input. 


Figure  4.  Experimental  setup.  A  commerce,  ultrasound  probe  is  locared  a,  die  top  of  the  water 
tank  and  the  MR  source  and  detector  fibers  are  deployed  around  the  tank  and  configured  m 


transmission  geometry. 


Figure  5.  Configurations  of  NIR  sources  and  detectors  used  in  the  reported  expenments 

(a)  Frequency  domain  experiments,  (b)  DC  experiments,  (c)  Top  view  of  target  rotation 
scheme  used  in  frequency  domain  measurements. 


Figure  6.  Comparison  of  reconstructed  absorption  coefficient  maps  of  a  high-contrast  target 
located  closer  to  die  boundary,  (a)  Reconstructed  p,distribution  using  NIR  data  only. 

(b)  Reconstructed  p,  distribution  with  ultrasound  localization,  <0  B-scan  ultiasonnd  image,  (d) 
Co-registered  C-scan  ultrasound  image,  (e)  Comparison  of  p,  distributions  along  x-axis,  pointed 
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by  two  vertical  arrows  in  part  (a),  (f)  Comparison  of  ^distributions  along  y-axis,  pointed  by  two 
horizontal  arrows  in  part  (a). 

Figure  7.  Comparison  of  reconstructed  absorption  maps  of  a  high-contrast  target  located  at 
the  middle  of  the  turbid  medium. 

(a)  Reconstructed  14  distribution  using  NIR  data  only,  (b)  Co-registered  C-scan 
ultrasound  image,  (c)  Reconstructed  |4  distribution  with  ultrasound  localization. 

Figure  8.  Comparison  of  the  reconstructed  absorption  maps  of  a  low-contrast  target  located 
closer  to  the  boundary,  (a)  Reconstructed  )4  distribution  using  NIR  data  only,  (b) 

Reconstructed  14 distribution  with  ultrasound  localization. 

Figure  9.  Comparison  of  reconstructed  absorption  coefficient  maps  of  a  high-contrast 
target  located  closer  to  the  boundary.  DC  system  of  configuration  Fig.5  (b)  was  used  for 

i 

experiments,  (a)  Reconstructed  )4  distribution  using  NIR  data  only,  (b)  Co-registered  C- 
scan  ultrasound  image,  (c)  Reconstructed  (4  distribution  with  ultrasound  localization,  (e) 
Comparison  of  |4  distribution  along  x-axis  as  pointed  by  vertical  arrows  shown  in  (c) .  (f) 
Comparison  of  J4  distribution  along  y-axis  as  pointed  by  horizontal  arrows  shown  in  (c). 
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Figure  10.  Comparison  of  reconstructed  absorption  maps  of  a  low-contrast  target  located 
closer  to  the  boundary.  DC  system  experiments,  (a)  Reconstructed  )J.a  distribution  using 

NIR  data  only,  (b)  Reconstructed  Jia  distribution  with  ultrasound  localization. 
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Table  I.  Experimental  results  obtained  from  the  frequency  domain  system 


Target  location 

NIR  only 

(Mean  fia  (cm-1) 

NIR  +  US 

(Mean  jla  (cm-1) 

True  value  | 
(Mean  (cm-1) 

Closer  to  the  boundary  (high  contrast) 

0.048(24%) 

0.149(75%) 

CM 

o 

(low  contrast) 

- 

0.078(78%) 

0.1 

At  the  middle  (high  contrast) 

0.052(26%) 

0.134(67%) 

0.2 

(low  contrast) 

-- 

0.072(72%) 

0.1 
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Table  II.  Experimental  results  obtained  from  the  DCsystem 


Target  location 

NIR  only 

(Mean  fla  (cm-i)) 

NIR  +  US 

(Mean//fl  (cm-i)) 

True  value 

(Mean  fla  (cm-i)) 

Closer  to  the  boundary  (high  contrast) 

0.0447(22%) 

0.153(77%) 

0.2 

(low  contrast) 

- 

0.087(87%) 

0.1 

At  the  middle  (high  contrast) 

0.040(20%) 

0.151(76%) 

0.2 

(low  contrast) 

- 

0.085(85%) 

0.1 
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Abstract 


In  this  letter,  a  novel  two-step  reconstruction  scheme  using  a  combined  near  infrared  and 
ultrasound  technique  and  its  utility  in  imaging  distributions  of  optical  absorption  and  hemoglobin 
concentration  of  breast  lesions  is  demonstrated.  In  the  first-step  image  reconstruction,  the  entire 
tissue  volume  is  segmented  based  on  initial  co-registered  ultrasound  measurements,  into  lesion 
region  and  background  region.  Reconstruction  is  performed  by  using  a  finer  grid  for  lesion 
region  and  a  coarse  grid  for  the  background  tissue.  As  results,  the  total  number  of  voxels  with 
unknown  absorption  can  be  maintained  on  the  same  order  of  total  measurements  and  the  matrix 
with  unknown  total  absorption  distribution  is  appropriately  scaled  for  inversion.  In  the  second- 
step,  image  reconstruction  is  refined  by  optimizing  lesion  parameters  measured  from  ultrasound 
images.  It  is  shown  in  this  letter  that  detailed  distributions  of  wavelength-dependent  absorption 
and  hemoglobin  concentration  of  breast  carcinoma  can  be  obtained  with  the  new  reconstruction 
scheme. 

OCIS  codes:  170.0170, 170.3010, 170.5270, 170.7170, 170.3830 


Tumor  blood  volume  and  micro-vascular  density  are  parameters  anatomically  and  functionally 
associated  with  tumor  angiogenesis.  If  a  single  wavelength  is  used,  optical  absorption  related  to 
angiogenesis  and  other  normal  blood  vessels  can  be  measured.  If  two  or  more  optical 
wavelengths  are  used,  oxy-hemoglobin  and  deoxy-hemoglobin  concentrations  can  be  measured 
simultaneously  and  total  hemoglobin  concentration  can  be  quantitatively  estimated.  During  the 
last  decade,  modeling  of  the  light  propagation  in  the  near  infrared  (NIR)  region,  combined  with 
the  advancements  of  light  source  and  detectors,  has  improved  the  diffused  light  measurements 
and  made  possible  the  application  of  tomographic  techniques  for  characterizing  and  imaging 
tumor  angiogenesis.  13  However,  the  NIR  technique  has  not  been  widely  used  in  clinics  and  the 
fundamental  problem  remains  the  intense  light  scattering.  As  a  result,  diffusive  light  probes  a 
widespread  region  instead  of  providing  information  along  a  straight  line,  and  tomographic  image 
reconstruction  is,  in  general,  underdetermined  and  ill-posed.  We  have  demonstrated  a  combined 
imaging  technique  using  a  priori  lesion  structure  information  provided  by  co-registered 
ultrasound  images  to  assist  NIR  imaging  reconstruction  in  phantom  studies.4'5  As  a  result,  the 
NER  image  reconstruction  is  well  defined  and  less  sensitive  to  noise.  We  have  also  demonstrated 
the  feasibility  of  using  NIR  imaging  as  adjunct  to  ultrasound  for  breast  cancer  diagnosis.6  In  this 
letter  we  report  our  novel  two-step  image  reconstruction  scheme  using  the  combined  approach 
and  demonstrate  its  utility  in  imaging  tumor  absorption  and  hemoglobin  distributions.  It  is 
shown  that  the  detailed  heterogeneous  distributions  of  wavelength-dependent  optical  absorption 
and  hemoglobin  concentration  of  a  breast  carcinoma  can  be  obtained.  To  the  best  of  our 
knowledge,  such  detailed  distributions  have  not  been  achieved  by  NIR  only  imaging  reported  in 


the  literature. 


A  picture  of  our  combined  hand-held  probe  used  in  clinical  studies  is  shown  in  Fig.  1(a),  and  the 
probe  dimensions  and  optical  sensor  distributions  are  shown  in  Fig.  1(b).  The  combined  probe 
consists  of  a  commercial  ultrasound  one-dimensional  array  located  at  the  center  of  the  probe  and 
optical  source  and  detector  fibers  distributed  at  the  periphery  and  connected  to  a  NIR  imager. 
The  NIR  imager  consists  of  12  dual-wavelength  source  channels  and  8  parallel  receiving 
channels.5  On  the  transmission  part,  12  pairs  of  dual  wavelength  (780nm  and  830nm)  laser 
diodes  are  amplitude  modulated  at  140  MHz.  On  the  reception  part,  8  photo  multiplier  tubes 
(PMTs)  detect  diffusely  reflected  light  from  the  tissue.  The  sources  and  detectors  are  coupled  to 
the  probe  through  optical  fibers.  Both  amplitude  and  phase  at  each  source-detector  pair  are 
obtained  and  the  resulting  total  number  of  measurements  is  12x8x2=  192.  The  combined 
probe  is  made  of  a  black  plastic  plate  10  cm  in  diameter,  therefore,  a  semi-infinite  boundary 
condition  can  be  used  for  NIR  measurement  geometry.  The  measured  amplitude  and  phase  after 
calibration  follow  simple  equations  as 

jlog(PapAap)  =  — kiPap  +  C, 

{^aP  =  ^rPajJ  +  ^2 

A  A 

where  amplitude  Aap  and  phase  are  associated  with  source  channel  a  and  detector  channel 
(3,  pap  is  the  corresponding  source-detector  separation,  k  is  the  wavenumber  given  by  k  = 

^  p~~ ^ ^  D  is  die  diffusion  coefficient  given  by  D=  l/(3fi’ ),  ©  and  v  are 

modulation  frequency  and  speed  of  light,  and  Cl  and  C2  are  constants.  By  fitting  the  amplitude 
and  phase  measured  from  the  normal  side  of  the  breast  we  can  obtain  kr  and  k;  and  calculate 
background  as 

pa  =  -D(kj  -  k? ) ,  j l\  =1/(3  D )  with  D = ©/(2  vkrk; ) 


(2) 


ARE 

MISSING 

IN 

ORIGINAL 

DOCUMENT 


In  our  two-step  image  reconstruction,  we  first  segment  tissue  volume  into  two  regions  L  and  B, 
which  contains  a  lesion  as  measured  from  co-registered  ultrasound  images  and  background 

tissue,  respectively.  We  use  Bom  approximation  to  relate  the  scattered  field  Usc(rsi,rdi,co) 
measured  at  the  source-detector  pair  i  to  absorption  variations  Apa(r’)  in  each  volume  element  of 
two  regions  within  the  sample 


UscUsi’^iii’®)  j-j 


jG(r\rdi)Uinc  (r\  rsi  )Apa  (r ’)d3r’+ J  G(r rdi  )Uinc  (r  ’ ,  rsi  )A(ia  (r ’)d3r  ’ 


(3) 


where  Uinc(r’,rsi,co)  and  G(r’,rdi,CD)  are  incident  wave  and  green  function  of  semi-infinite 
geometry,  respectively.  rsi  and  rdi  are  source  and  detector  positions.  We  then  discretize  the 

lesion  volume  and  background  volume  with  different  voxel  size  (a  finer  grid  for  lesion  volume 
and  a  coarse  grid  for  background).  The  scattered  field  can  then  be  approximated  as 


U*  Osi,rdi 


£G(rvj,rdi  )Uinc(rVj,rsi)|  Apa  (r’)d3r’+£G(rvk  ,rdi  )Uioc(rvk,rsi  )J  Ajia  (r’)dV 

Lj  j  Bk  k 


(4) 


where  Gfr^r^U^fr^rJ  and  GCr^r^U^r*^)  are  weight  values  at  the  center  rvjof 
voxel  j  and  center  rvk  of  voxel  k  in  lesion  volume  L  and  background  volume  B,  respectively. 
The  matrix  form  of  equation  (4)  is  given  as 

[U«,U=[wl.wbLn[ml.Mb]t  (5) 


where  W^  —  [  ^  G(rVj ,  rdj  )U  -nc  (rVj ,  rS[  and  ^  G(rvk ,  rd;  )Umc  (rvk ,  rSj  )]f^XNg  ^r£ 

weight  matrixes  for  lesion  volume  and  background  volume,  respectively; 

[Ml]  =  [ J A|ia(r’)d3r\ . J AM-a(r’)d3r1  and  [MB]  =  [|Aiia(r’)dV . |Am(r’)dV]  are 

il  nl  ib  nb 


total  absorption  distributions  of  lesion  volume  and  background  volume ,  respectively. 


Instead  of  reconstructing  Apa  distribution  directly  as  the  standard  Bom  approximation,  we 
reconstruct  total  absorption  distribution  M  and  then  divide  the  total  by  different  voxel  sizes  of 
lesion  and  background  tissue  to  obtain  Apa  distribution.  By  choosing  a  finer  grid  for  lesion  and 
a  coarse  grid  for  background  tissue,  we  can  maintain  the  total  number  of  voxels  with  unknown 
absorption  on  the  same  scale  of  the  total  measurements.  As  a  result,  the  inverse  problem  is  less 
underdetermined  and  ill-posed.  In  addition,  since  the  lesion  absorption  coefficient  is  higher  than 
that  of  background  tissue,  in  general,  the  total  absorption  of  the  lesion  over  a  smaller  voxel  is  on 
the  same  scale  of  total  absorption  of  the  background  over  a  bigger  voxel,  therefore  the  matrix 
[Ml,Mb]  is  appropriately  scaled  for  inversion.  The  reconstruction  is  formulated  as  least  square 

problem  of  minimizing  object  function  min  ||usd  —  Wlm£  -  WBMB  .  The  unknown 
distribution  M  can  be  iteratively  calculated  using  conjugate  gradient  search  method. 

The  lesion  location  and  volume  obtained  from  co-registered  ultrasound  images  are  estimated  as 
follows.  Since  the  commercial  1-D  ultrasound  probe  we  used  acquires  2-D  ultrasound  images  in 
y-z  plane  (z  is  the  propagation  direction)  and  the  2-D  NIR  probe  provides  3-D  images,  the  co¬ 
registration  is  limited  to  an  interception  plane.  However,  if  we  approximate  a  lesion  as  an 
ellipsoid,  we  are  able  to  estimate  its  radii  from  two  orthogonal  ultrasound  images  and  therefore 
to  obtain  its  volume.  The  3-D  lesion  center  is  approximated  by  first  obtaining  a  2-D  lesion 
center  from  a  2-D  ultrasound  image  and  then  perturbing  the  2-D  center  in  the  other  orthogonal 
spatial  direction  in  the  second  step  reconstruction  discussed  below  in  an  attempt  to  find  out  the 
best  estimate  of  the  center  coordinate  in  that  direction.  Another  source  of  error  in  estimating  the 


lesion  geometry  is  the  inaccuracy  in  radius  measurements  because  lesion  boundaries  may  not  be 
well  defined  in  ultrasound  images.  In  addition,  the  target  volumes  seen  by  different  modalities 
may  be  different  because  of  different  contrast  mechanisms.  In  the  second  step,  we  refine  image 
reconstruction  by  perturbing  the  center  c0  =  (x0,y0,z0)  and  then  the  radii  r0  =  (rx,ry,zy)  and 

choosing  the  optimal  set  of  parameters  ( copI ,  ropt).  The  final  refined  reconstruction  is 
formulated  as 


min||uSd  -  WL(c 

Opt  >  *opt  )ML  (C0pt  >  ropt  )  (C0pt  >  ^opt  L  (^opt  ’  *opt  )| 


(6) 


Conjugate  graduate  method  is  used  to  search  for  global  minimum.  Since  the  total  absorption 
matrix  M  is  appropriately  scaled,  the  iterative  search  converges  quickly  in  two  to  three  iterations. 
Two  criteria  were  used  for  choosing  (copt,  ropt).  First,  the  initial  radii  measured  from  co¬ 
registered  ultrasound  were  fixed  and  optimal  center  copt  was  chosen  which  produced  maximum 

total  absorption  energy  inside  the  lesion.  Second,  the  optimal  center  was  used,  and  radii  were 
increased  or  decreased  to  obtain  optimal  radii,  which  produced  maximum  total  absorption  energy 
inside  the  lesion. 


Clinical  studies  were  performed  at  the  Health  Center  of  the  University  of  Connecticut  and  the 
human  subject  protocol  was  approved  by  Health  Center  IRB  committee.  Patients  with  palpable 
and  non-palpable  masses  that  were  visible  on  clinical  US  were  used  as  subjects.  These  subjects 
were  scanned  with  the  combined  probe,  and  ultrasound  images  and  optical  measurements  were 
acquired  at  multiple  locations  including  the  lesion  region  scanned  at  two  orthogonal  positions,  a 


normal  region  of  the  same  breast  scanned  at  two  orthogonal  positions  (if .the  breast  is  large)  and  a 
normal  region  of  the  contralateral  breast  scanned  at  two  orthogonal  positions. 


An  example  is  given  in  this  letter  to  demonstrate  the  use  of  our  reconstruction  scheme.  Figure  2 
(a)  shows  a  gray  scale  ultrasound  image  of  a  palpable  lump  of  a  44-year-old  woman.  The  lesion 
was  located  at  6  to  8  o’clock  position  of  the  left  breast  at  approximately  1.5  cm  depth. 
Ultrasound  showed  an  irregular  poorly  defined  hypoechoic  mass  and  the  lesion  was  considered 
as  highly  suspicious  for  malignancy.  An  ultrasound  guided  core  needle  biopsy  was 
recommended.  Biopsy  results  revealed  that  the  lesion  was  a  high  grade  in-situ  ductal  carcinoma 
with  necrosis. 


Multiple  optical  measurements  at  two  orthogonal  positions  were  simultaneously  made  with 
ultrasound  images  at  the  lesion  location  as  well  as  at  approximately  the  same  location  of  the 
contralateral  normal  breast.  The  fitted  average  tissue  background  measured  at  normal  side  of  the 

breast  at  both  wavelengths  were  |ia8°=  0.03  cm'1,  |ia  =  0.05  cm'1,  |is  =9.22  cm'1, 
|l’  =7.61  cm'1.  The  perturbations  for  both  wavelengths  used  to  calculate  absorption  maps 


were  normalized  as 


where 

UN(rsi,rdi,©) 


UL(rsi,rdj,co)  and  UN(rsi,rdi,co)  were  measurements  obtained  at  lesion  region  and  contralateral 
normal  region,  and  UB(rsi,rdi,G>)  was  calculated  incident  field  using  fitted  background  (J.aand 

D=l/3p’  .  This  procedure  cancels  unknown  system  gains  associated  with  different  sources  and 
detectors  as  well  as  electronic  channels.  The  initial  estimate  of  lesion  center  and  diameters  from 


ultrasound  images  were  (0, 0.39  cm,  1.7  cm)  and  3.44  cm  x  4.38cm  x  1.76  cm,  respectively.  A 
finer  grid  of  0.5  x  0.5  x  0.5  (cm3)  and  a  coarse  grid  of  1.5  x  1.5  x  1.0  (cm3)  were  chosen  for  the 
lesion  and  background  tissue,  respectively.  The  total  reconstruction  volume  was  chosen  to  be  9 
x  9  x  4  cm3  and  the  total  voxels  with  unknown  optical  absorption  was  190,  which  was  on  the 
same  order  of  192  total  measurements.  The  image  reconstruction  was  performed  using  the  NIR 
data  simultaneously  acquired  with  the  ultrasound  image  shown  in  Fig.2  (a).  The  second  step 
refined  reconstruction  revealed  optimal  lesion  centers  at  approximately  (-1.1  cm,  0.3  cm,  1.7  cm) 
for  780  nm  and  (-0.9  cm,  -0.7  cm,  1.7  cm)  for  830  nm,  respectively,  and  optimal  diameters  of 
4.28  cm  x  5.18  cm  x  1.96  cm.  The  detailed  absorption  maps  with  high  absorption  non-uniformly 
distributed  around  the  lesion  boundaries  at  both  wavelengths  are  shown  in  Fig.  2  (b)  and  (c).  By 
assuming  that  the  major  chromophores  are  oxygenated  (oxyHb)  and  deoxygenated  (deoxyHb) 
hemoglobin  molecules  in  the  wavelength  range  studied,  we  can  estimate  the  distribution  of  total 
hemoglobin  concentration  as  shown  in  Fig.2  (d).  The  extinction  coefficients  used  for 
calculating  oxyHb  and  deoxyHb  concentrations  were  =2.544,  £^°02=  1.695,  e^b  =1-797-, 

eHb°02=2-419  obtained  in  Ref.  7  in  a  natural  logarithm  scale  with  units  of  inverse  millimoles 

times  inverse  centimeters.  The  measured  average  cancer  and  background  total  hemoglobin 
concentrations  were  1 16.9  |i  moles  and  47.6  n  moles,  respectively.  These  values  agree  with  3-D 
average  measurements  reported  in  Ref.  8. 

It  is  interesting  to  note  that  the  absorption  distributions  at  both  wavelengths  as  well  as  total 
hemoglobin  concentration  were  distributed  heterogeneously  at  the  cancer  periphery.  To  the  best 
of  our  knowledge,  such  fine  distributions  have  not  been  reported  by  using  NIR  only 
reconstruction  techniques.  However,  this  finding  agrees  with  the  published  literature  showing 


that  breast  cancers  have  higher  blood  volumes  than  non -malignant  tissue  due  to  angiogenesis, 
especially  at  the  cancer  periphery.9  In  addition,  the  carcinoma  reported  here  had  a  necrotic  core 
which  could  lead  to  low  absorptions  observed  at  both  wavelengths  in  the  center  region.  Ideally, 
pathologic  report  with  blood  vessel  density  count  could  be  used  to  correlate  the  NIR  findings 
reported  here.  Future  plans  are  to  obtain  microvessel  counts  in  the  region  of  the  tumor  as  well  as 
surrounding,  non-malignant  breast  tissue. 


In  summary,  we  have  reported  a  two-step  image  reconstruction  scheme  and  we  demonstrate  that 
aided  with  co-registered  ultrasound,  detailed  optical  absorption  and  total  hemoglobin 
distributions  of  the  breast  carcinoma  can  be  obtained.  These  distributions  reveal  tumor 
angeogenesis  and  provide  valuable  information  for  breast  cancer  diagnosis  and  treatment. 
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Fig.  1.  (a)  Picture  of  a  hand-held  combined  probe  and  a  frequency  domain  NIR  imager,  (b) 

Sensor  distribution  of  the  combined  probe.  The  diameter  of  the  combined  probe  is  10  cm. 
Smaller  circles  are  optical  source  fibers  and  big  circles  are  detector  fibers.  A  commercial 
ultrasound  probe  is  located  at  the  center  and  its  dimensions  are  5.6  cm  by  1  cm. 


Fig.  2.  (a)  A  gray  scale  ultrasound  image  of  a  palpable  lump  of  a  44-year-old  woman. 
Ultrasound  showed  an  irregular  poorly  defined  hypoechoic  mass  and  the  lesion  was  considered 
as  highly  suspicious  for  malignancy.  Reconstructed  optical  absorption  maps  at  780  nm  (b)  and 
830  nm  (c).  Vertical  color  bars  are  absorption  coefficient  in  the  unit  of  cm'1,  (d)  is  the  total 
hemoglobin  concentration  map.  The  vertical  color  bars  are  p.  moles.  The  NIR  data  were 
simultaneously  acquired  with  ultrasound  image  shown  in  Fig.2  (a).  Each  image  consists  of  7 
slices  obtained  in  0.5  cm  spacing  from  0.5  cm  to  3.5  cm  in  depth.  The  vertical  and  horizontal 
axes  correspond  to  x  and  y  dimensions  of  9  cm  by  9  cm. 
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as  highly  suspicious  for  malignancy.  Reconstructed  optical  absorption  maps  at  780  nm  (b)  and 
830  nm  (c).  Vertical  color  bars  are  absorption  coefficient  in  the  unit  of  cm'1,  (d)  is  the  total 
hemoglobin  concentration  map.  The  vertical  color  bars  are  p  moles.  The  NIR  data  were 
simultaneously  acquired  with  ultrasound  image  shown  in  Fig.2  (a).  Each  image  consists  of  7 
slices  obtained  in  0.5  cm  spacing  from  0.5  cm  to  3.5  cm  in  depth.  The  vertical  and  horizontal 
axes  correspond  to  x  and  y  dimensions  of  9  cm  by  9  cm. 
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We  propose  a  novel  noniterative  near-infrared  diffusive  image  reconstruction  method  that  uses  minimal 
a  priori  co-registered  ultrasound  information.  Small  absorbing  targets  embedded  in  a  homogeneous 
background  are  described  approximately  in  terms  of  their  monopole,  dipole,  and  quadrupole  moments.  With 
an  approximate  estimation  of  the  center  locations  of  these  absorbers  from  ultrasound  images,  we  show  in 
simulations  that  the  reconstruction  accuracy  of  the  absorption  coefficient  exceeds  80%  if  the  noise  level  is 
less  than  0.2%.  We  also  demonstrate  experimentally  that  the  accuracy  can  be  improved  by  use  of  additional 
ultrasound  volume  information  even  for  a  noise  level  as  high  as  1.5%.  ©  2002  Optical  Society  of  America 
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Functional  imaging  with  near-infrared  (NIR)  light 
has  found  potential  applications  in  many  areas, 
such  as  breast  and  brain  lesion  detection  and  diag¬ 
nosis.1,2  Recently,  a  combination  of  NIR  imaging 
with  other  imaging  modalities,  such  as  ultrasound  or 
magnetic  resonance  imaging,  has  shown  promising 
results3"5  in  providing  complementary  contrasts  and 
overcoming  NIR  reconstruction  problems  related  to 
intensive  light  scattering. 

In  this  Letter  we  introduce  a  novel  reconstruction 
algorithm  for  NIR  diffusive  imaging  that  uses  approxi¬ 
mate  target  center  locations  estimated  from  co¬ 
registered  ultrasound.  Our  new  method  is  based 
on  estimation  of  major  characteristics  of  isolated 
small  absorbers.  These  characteristics  are  monopole, 
dipole,  and  quadrupole  moments.  Higher-order  mo¬ 
ments  have  negligible  effects  in  characterizing  the 
absorbers  and  are  ignored.  It  is  shown  in  this  Letter 
that  measurements  of  diffusive  photon  density  waves 
cannot  readily  achieve  a  signal-to-noise  ratio  that  is 
good  enough  for  reconstructing  the  detailed  shape 
information  of  targets.  However,  the  quadrupole 
moment  can  provide  an  approximate  extension  of  the 
target,  which  is  necessary  for  distributing  the  integral 
absorption  (the  monopole)  to  an  appropriate  target  re¬ 
gion.  In  situations  in  which  the  measurement  system 
suffers  from  unexpected  noise,  accurate  estimation  of 
the  quadrupole  moment  might  not  be  possible.  The 
target  volume  estimated  from  co-registered  ultrasound 
images  can  be  used  for  distributing  the  integral  ab¬ 
sorption  inside  the  target  volume.  We  have  conducted 
a  series  of  simulations  at  different  noise  levels  to 
evaluate  the  performance  of  this  new  method.  We 
have  also  conducted  experiments  using  our  combined 
NIR-ultrasound  imager5  to  test  the  algorithm. 

If  it  is  assumed  that  there  are  only  a  few  small  inho¬ 
mogeneous  targets  embedded  in  a  homogeneous  back¬ 
ground,  the  distribution  of  the  absorption  coefficient 
can  be  expressed  as 

lia(r)  =  fia°  +  A lia(r) .  (1) 

The  background  absorption  coefficient,  /aq°,  can  be  as 
small  as  0.02  cm-1  for  normal  breast  tissues,  whereas 
A fia  can  be  well  beyond  0,1  cm"1  for  tumors  because 
of  the  blood.  In  the  Born  approximation,  the  photon 
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density  wave  is  decomposed  into  an  incident  wave  and 
a  scattering  wave: 

<P(<o,r)  =  $inc(*M’)  +  (2) 

The  incident  wave  is  the  solution  of  a  homogeneous  dif¬ 
fusion  equation,  and  the  scattering  wave  can  be  given 
approximately  by  the  linear  perturbation  theory: 

$9Ct(w,r)=  -  <i>inc(^y)G(r,r')rfV.  (3) 

In  this  equation,  G(r,r')  is  the  Green’s  function, 
and  D  =  1/3 is  inversely  proportional  to  the  re¬ 
duced  scattering  coefficient,  jig1.  In  the  context  of 
early  stage  breast  cancer  detection,  we  can  further 
assume  that  those  heterogeneities  are  confined  to  a 
few  insolated  regions.  Then,  the  scattering  wave  can 
be  expanded  around  each  of  N  target  centers  as 

1  N 

*«*(«, r)  =  -  jr  X  [M-Wtro")  +  D"  •  VW(r „*) 

+  Q1'  *  VVW(r0*/)/2  4-  0(a3)].  (4) 

Here  rov  is  the  center  of  the  yth  target,  and  W(ro")  — 
^incC^ro^Gforo1')  is  the  weight  function.  For  the 
yth  target,  M*,  D1',  and  CT  are  its  monopole,  dipole,  and 
quadrupole  moments,  respectively.  Terms  beyond  the 
second  order  of  target  dimension  a  have  negligible  ef¬ 
fects  in  characterizing  the  absorbers  and  are  neglected 
in  our  model.  For  the  pth  target,  its  monopole,  dipole, 
and  quadrupole  moments  are  scalar,  vector,  and  sec¬ 
ond-order  tensor,  respectively,  and  they  are  given  as 

M"  =  f  Afia(r')d3r' ,  D"  =  f  Ava(r')r'd: V , 

Q"  =  f  A/UrVr'dV.  (5) 

Jv„ 

The  integrations  are  over  a  small  isolated  region  Vv. 

A  multiple-source,  multiple-detector  configuration 
is  typical  for  frequency-domain  diffusive  imaging 
systems.  We  are  using  12  sources  (780  nm)  that  are 
amplitude  modulated  at  140  MHz  and  eight  detectors 
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•in  our  simulations  as  well  as  in  the  experiments. 

.  These  sources  and  detectors  are  deployed  on  an  ab¬ 
sorbing  plane,  which  simplifies  boundary  conditions. 
The  reflection  mode  with  a  semi-infinite  geometry 
is  used  for  both  simulation  and  experiment.  The 
measured  scattering  wave  is  related  to  the  target 
moments  by  the  following  equation: 

M1  D1  Q1  ■■■  Qn]t.  (6) 

An  element  of  SI  is  either  the  weight  function  for  one 
target  or  its  derivatives  up  to  the  second  order. 

The  inverse  problem  is  to  retrieve  the  character¬ 
istics  of  those  embedded  absorbers  from  measured 
photon  density  waves  on  the  surface.  Based  on  our 
simplified  forward  model,  it  is  necessary  to  have 
certain  information  about  absorbers,  such  as  the  num¬ 
ber  of  targets  and  their  center  positions,  in  advance. 
The  ultrasound  system  of  our  combined  imager  can 
provide  such  information.  Once  the  initial  center 
locations  have  been  specified,  the  weight  function, 
together  with  its  first  and  second  derivatives,  can  be 
calculated  immediately.  Since  the  inverse  problem 
has  only  a  few  unknowns,  reconstruction  of  target 
moments  is  overdetermined  and  least-squares  solu¬ 
tions  can  be  adopted.  For  each  target,  the  monopole 
moment  represents  the  integral  absorption,  the  dipole 
moments  result  in  correction  of  the  center  position, 
and  the  quadrupole  moments  lead  to  estimation  of 
the  target  volume,  in  which  the  integral  absorption 
can  be  redistributed.  The  target  volume  is  estimated 
by  simple  matching  of  the  quadrupole  moments  with 
those  of  an  ellipsoid,  which  has  the  freedom  to  rotate 
about  the  center  in  any  direction  to  any  angle. 

Figure  1  shows  a  schematic  of  the  combined  probe 
for  experiments,  which  is  a  circular  plate  made  of  rigid 
absorbing  material.  One  or  two  cuboid  absorbers  were 
placed  2.5  cm  deep  in  a  homogeneous  medium.  For 
phantom  experiments,  raw  data  were  acquired  with 
our  combined  imaging  system.  The  ultrasound  array 
was  translated  in  the  x  direction  to  yield  target  vol¬ 
ume  estimates.  We  used  0.6%  Intralipid  as  a  homoge¬ 
neous  background,  which  yielded  a  reduced  scattering 
coefficient  of  6  cm-1  and  an  absorption  coefficient  of 
0.015  cm”1.  The  embedded  targets  were  as  scatter¬ 
ing  as  the  background  but  were  more  absorbing.  For 
simulations,  the  simulated  measurement  data  were  the 
sum  of  forward  model  prediction  and  random  noises 
ranging  from  0.1%  to  1%  (with  respect  to  the  ampli¬ 
tudes  of  the  incident  waves).  The  optical  properties 
of  homogeneous  media  and  targets  are  similar  to  those 
used  in  phantom  experiments. 

Shown  in  Fig.  2  are  simulation  results  for  a  single 
target.  The  target  was  a  cuboid  of  (1.2,  0.8,  1)  cm 
in  the  x ,  y,  and  2  directions,  respectively,  and  was 
centered  at  (0.2,  0,  2.5)  cm.  The  volume  of  the  target 
was  0.96  cm3.  The  absorption  coefficient  above  the 
background  was  0.25  cm-1.  Here  we  provide  the 
results  for  one  dipole  moment  (Dx)  and  one  quadruple 
moment  ( Qxx ),  but  similar  results  were  obtained  for 
other  components.  The  initial  target  position  that 
was  input  into  the  reconstruction  algorithm  was  (0, 


0,  2.5)  cm,  which  deviated  slightly  from  the  true 
position.  This  initial  location  error  resulted  in  a 
negligible  effect  on  the  reconstructed  monopole  values, 
as  shown  in  Fig.  2(a).  The  offset  of  center  position 
can  be  corrected  by  the  estimated  dipole  moment, 
Dx  [see  Fig.  2(c)],  which  is  slightly  lower  than  the 
true  value.  However,  the  inaccuracy  of  the  initial 
center  position  obviously  affects  the  reconstructed 
quadrupole  moment.  The  circles  in  Fig.  2(e)  repre¬ 
sent  values  greater  than  the  true  values.  After  we 
correct  the  center  position  according  to  the  estimated 
dipole  moments,  the  reconstructed  values  of  Q*x,  plot¬ 
ted  as  asterisks  in  the  same  figure,  are  much  closer  to 
the  true  values.  Although  the  mean  values  of  recon¬ 
structed  moments  are  essentially  independent  of  noise 


Fig.  1.  (a)  Bottom  view  of  the  combined  probe.  The 
central  rectangular  slot  gives  the  ultrasound  array  access 
to  tissues  underneath  the  probe.  The  circular  holes 
are  used  to  hold  optical  fibers.  The  small  holes  are  for 
light  sources,  and  the  larger  ones  are  for  detectors.  The 
diameter  of  the  probe  is  10  cm.  (b)  Side  view  of  the  probe. 


Fig.  2.  Simulation  results  for  a  single  target,  (a),  (c),  and 
(e)  mean  values  of  monopole,  dipole,  and  quadrupole  mo¬ 
ments,  respectively,  (b),  (d),  and  (f)  standard  deviations 
of  the  corresponding  quantities.  The  solid  lines  in  (a),  (c), 
and  (e)  are  true  values. 
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Table  1.  Reconstructed  Absorption  Coefficients  from  Phantom  Experiments0 


- - - - - - 

Combination  1 

Combination 

2 

Value  (cm"1) 

Target  1 

Target  2 

Target  1 

Target  2 

True 

Reconstructed 

0.065 

0.064  ±  0.028 

0.065 

0.064  ±  0.027 

0.102 

0.106  ±  0.026 

0.065 

0.064  ±  0.017 

°The  average  values  and  standard  deviations  are  estimated  from  six  sets  of  experimental  data. 


level,  the  standard  deviations  of  all  moments  increase 
linearly  with  increasing  noise  level.  Higher-order 
moments  suffer  more  from  noise,  as  expected  intui¬ 
tively.  For  the  quadrupole  moments,  as  small  as  a 
0.2%  noise  level  could  lead  to  a  relative  error  of —20%. 
This  is  the  main  reason  that  we  did  not  attempt  to 
reconstruct  moments  beyond  the  second  order. 

To  test  the  capability  of  our  algorithm  to  charac¬ 
terize  multiple  targets,  we  performed  simulations  with 
two  absorbers  of  different  absorption  coefficients.  The 
absorbers  were  the  same  size  as  the  single  target.  One 
absorber  was  located  at  (-1.2, 0,2.5)  cm,  with  an  ab¬ 
sorption  coefficient  0.25  cm"1  beyond  the  background. 
The  other  absorber  was  located  at  (1.2, 0, 2.5)  cm,  with 
an  absorption  coefficient  0.1  cm-1  beyond  the  back¬ 
ground.  The  main  effect  of  adding  one  target  was  that 
the  relative  error  in  the  reconstructed  monopole  mo¬ 
ments  increased  to  —2.5  times.  Less  of  an  effect  on  the 
dipole  moments  was  observed,  although  an  even  more 
trivial  difference  for  the  quadrupole  could  be  seen.  It 
was  also  found  that  the  noise-induced  deviation  was 
independent  of  target  optical  properties. 

Our  current  NIR  imaging  system  has  a  noise  level 
near  1.5%.  According  to  the  simulation  results,  the 
system  cannot  readily  retrieve  the  accurate  quadrupole 
moments  from  measured  photon  density  waves.  Thus 
we  have  to  use  the  volume  information  obtained  from 
co-registered  ultrasound  images  in  target  characteri¬ 
zation.  We  used  two  geometrically  identical  targets 
that  were  approximately  1-cm3  cubes.  They  were  gel 
phantoms  made  from  0.6%  Intralipid  solution,  ink, 
and  ultrasound  scatterers.  The  boundaries  of  these 
targets  can  clearly  be  seen  in  ultrasound  images,  and 
target  centers  and  volumes  can  be  estimated  accu¬ 
rately.  For  example,  a  1.07-cm3  target  was  measured 
as  1.05  cm3  by  ultrasound.  For  co-registered  NIR 
imaging,  two  different  combinations  of  targets  were 
adopted.  The  first  was  a  pair  of  equally  absorbing 
cubes  with  an  absorption  coefficient  of  0.065  cm'1. 
The  second  was  a  pair  of  absorbing  cubes  with  absorp¬ 
tion  coefficients  of  0.065  and  0.1  cm"1,  respectively. 
We  simply  distribute  the  estimated  monopole  evenly 
over  the  target  volume.  So  the  reconstructed  ab¬ 
sorption  coefficient  was  an  average  value  within  each 
target.  The  true  and  the  reconstructed  values  are 
compared  in  Table  1.  The  mean  values  and  standard 
deviations  were  based  on  six  sets  of  measurements. 

Since  target  monopole  moments  can  be  estimated 
with  sufficient  accuracy,  the  reconstructed  absorption 
coefficient  depends  more  on  the  estimation  of  target 
volume,  which  can  be  obtained  from  the  reconstructed 
quadrupole  moment  or  co-registered  ultrasound  im¬ 
ages.  It  is  possible  that  the  extent  of  the  optical 
volume  of  a  lesion  is  larger  than  the  volume  esti¬ 


mated  from  the  acoustic  images.  This  is  one  of  our 
motivations  for  estimating  the  optical  volume  from 
the  quadrupole  moments.  Further  improvements 
are  needed  for  an  accurate  estimate  of  quadrupole 
moments.  In  the  moment-based  method  we  assume 
that  the  lesions  are  isolated  and  are  embedded  in  a 
homogeneous  background.  This  assumption  is  quite 
true  for  more-homogeneous  fatty  breasts  and  may  not 
hold  for  dense  breasts  that  consist  of  both  glandular 
tissue  and  fat.  In  the  latter  case,  we  could  segment 
ultrasound  images,  identify  tissue  types,  and  estimate 
background  optical  properties  of  different  tissues  in 
the  reconstruction  as  well.  Both  acoustic  and  optical 
contrasts  exist  in  tumors  but  the  sensitivities  of  these 
two  modalities  may  be  different.  Therefore,  the  cor¬ 
relation  between  acoustic  and  optical  heterogeneities 
remains  unknown  at  the  current  stage  and  will  be 
determined  in  our  future  clinical  studies. 

-  In  comparison  with  current  sophisticated  NIR  image 
reconstruction  algorithms,  our  method  is  very  simple 
in  terms  of  computation.  Once  an  ultrasound  image 
has  been  taken  and  targets  have  been  identified,  it 
is  possible  to  monitor  the  local  change  of  absorption 
continuously  and  in  real  time.  Averaging  over  a  large 
number  of  samples  will  help  suppress  the  effect  of  ran¬ 
dom  noise. 

In  this  Letter  an  algorithm  for  characterizing 
absorbers  has  been  demonstrated.  Theoretically,  a 
similar  algorithm  is  possible  for  characterizing  scat¬ 
tering  targets  as  well.  However,  since  the  weights 
of  monopole,  dipole,  and  quadrupole  moments  of 
scattering  targets  are  related  to  spatial  derivatives, 
which  are  two  order  higher  than  those  of  absorbers, 
the  quadrupole  moments  will  be  even  harder  to 
reconstruct. 
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Near  Infrared  Diffusive  Light 
Imaging  with  Ultrasound 
Localization 

Quing  Zhu,  Nan  Guang  Chen,  Puyun  Guo, 

Shikui  Yan,  and  Daqing  Piao 

Functional  imaging  with  near  infrared 
(NIR)  diffusive  light  has  found  poten¬ 
tial  applications  in  many  areas.  Although 
huge  efforts  have  been  made  to  bring  this 
technique  to  a  clinically  acceptable  stage, 
intense  light  scattering  remains  a  funda¬ 
mental  problem.  Subsequently,  diffusive 
light  continues  to  probe  a  widespread  re¬ 
gion  instead  of  providing  information 
along  a  straight  line.  Multiple  measure¬ 
ments  are  always  correlated  as  a  result  of 
overlapping  probed  regions.  Therefore, 
increasing  the  total  number  of  measure¬ 
ments  does  not  necessarily  provide  more 
independent  information  for  image  re¬ 
construction.  Image  reconstruction  is 
generally  affected  by  many  factors,  includ¬ 
ing  system  signal  to  noise  ratio,  probe  con¬ 
figuration,  and  regulation  schemes  used  in 
the  imaging  algorithms. 

We  have  constructed  a  novel  combined 
imager  suitable  for  simultaneous  NIR  dif¬ 
fusive  light  and  ultrasound  imaging  and 
co-registration. l'2  The  imager  consists  of  a 
two-dimensional  (2D)  ultrasound  array 
deployed  at  the  center  of  a  hand-held 
probe.  Twelve  dual  wavelength  laser  source 
fibers  and  eight  optical  detector  fibers  are 
deployed  at  the  periphery  of  the  probe. 
The  combined  imaging  uses  target  loca¬ 
tions  and  shapes  provided  by  co-registered 
ultrasound  images  to  localize  the  NIR  im¬ 
aging  reconstruction.  As  a  result,  the  NIR 
image  reconstruction  is  well  defined. 
Since  the  convergence  can  be  achieved 
within  a  small  number  of  iterations,  the 
NIR  image  reconstruction  is  less  sensitive 
to  noise. 
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In  experiments  where  co-registered  ul¬ 
trasound  provided  the  target  depth  infor¬ 
mation,  the  reconstructed  optical  absorp¬ 
tion  coefficient  accuracy  has  been  im¬ 
proved  by  15%  for  high  contrast  agents 
and  30%  for  low  contrast  agents3  (AQ:1). 
The  speed  of  reconstruction  has  been  im¬ 
proved  by  ten  times  on  average.  In  exper¬ 
iments  where  the  three-dimensional  (3D) 
target  distributions  were  provided  by  co¬ 
registered  ultrasound,  the  optical  image 
reconstruction  converges  very  fast  and 
needs  only  one  iteration  to  reconstruct  ac¬ 
curate  optical  absorption  coefficient.2 
This  result  is  very  encouraging  because 
there  is  no  known  robust  stopping  criteri¬ 
on  to  terminate  iterative  image  recon¬ 
struction  algorithms.  No  further  itera¬ 
tions  are  needed  with  the  a  priori  knowl¬ 
edge  from  co-registered  ultrasound. 

Figure  1  shows  an  example  of  ultra¬ 
sound  assisted  NIR  image  reconstruction.2 
Figure  1  (a)  shows  the  contours  of  ultra¬ 
sound  image  of  two  targets  of  1  cm  in  size. 
Figure  1  (b)  is  the  NIR  image  of  the  two 
targets.  The  target  pointed  by  the  wh]te_ 
arrow  in  part  (b)  is  displaced  from  its  true 
location  by  about  1  cm  and  the  recon¬ 
structed  absorption  coefficient  is  only 
50%  of  the  true  value.  With  the  ultra¬ 
sound  localization  shown  in  part  (a) ,  the 
NIR  image  shown  in  Fig.  1  (c)  provides 
correction  target  locations  and  the  recon¬ 
structed  absorption  coefficients  are  within 
20%  of  the  true  values.  By  adding  optical 
contrasts  to  ultrasonically  detected  lesions, 
the  combined  approach  has  great  promise 
in  improving  breast  cancer  diagnosis.4 
Optical  contrasts  will  significantly  im¬ 
prove  diagnostic  ultrasound  in  differenti¬ 
ating  benign  from  malignant  solid  lesions. 
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Figure  1.  (a)  Contours  of  the  ultrasound  image 
of  two  low  optical  contrast  targets  (absorption 
coefficients  are  0.1  cm-1)  embedded  2.5  cm  deep 
in  the  Intralipid  of  0.5%  concentration,  (b)  NIR  im¬ 
age  of  the  two  targets  without  ultrasound  localiza¬ 
tion,  One  target  (pointed  by  the  white  arrow)  is 
only  50%  of  the  true  value  and  is  displaced  from  its 
true  location  by  1  cm.  Two  image  artifacts  appear 
at  the  top  edge  of  the  image,  (c)  NIR  image  assist¬ 
ed  by  ultrasound  localization.  The  targets  appear 
at  correct  location  and  the  reconstructed  optical 
absorption  coefficients  are  within  20%  of  the  true 
value. 
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Simultaneous  near-infrared  diffusive  light  and 
ultrasound  imaging 


Nan  Guang  Chen,  Puyun  Guo,  Shikui  Yan,  Daqing  Piao,  and  Quing  Zhu 


We  have  constructed  a  near-real-time  combined  imager  suitable  for  simultaneous  ultrasound  and  near- 
infrared  diffusive  light  imaging  and  coregistration.  The  imager  consists  of  a  combined  hand-held  probe 
and  the  associated  electronics  for  data  acquisition.  A  two-dimensional  ultrasound  array  is  deployed  at 
the  center  of  the  combined  probe,  and  12  dual-wavelength  laser  source  fibers  (780  and  830  nm)  and  8 
optical  detector  fibers  are  deployed  at  the  periphery.  We  have  experimentally  evaluated  the  effects  of 
missing  optical  sources  in  the  middle  of  the  combined  probe  on  the  accuracy  of  the  reconstructed  optical 
absorption  coefficient  and  assessed  the  improvements  of  a  reconstructed  absorption  coefficient  with  the 
guidance  of  the  coregistered  ultrasound.  The  results  have  shown  that,  when  the  central  ultrasound 
array  area  is  in  the  neighborhood  of  2  cm  x  2  cm,  which  corresponds  to  the  size  of  most  commercial 
ultrasound  transducers,  the  optical  imaging  is  not  affected.  The  results  have  also  shown  that  the 
iterative  inversion  algorithm  converges  quickly  with  the  guidance  of  a  priori  three-dimensional  target 
distribution,  and  only  one  iteration  is  needed  to  reconstruct  an  accurate  optical  absorption  coefficient. 
©  2001  Optical  Society  of  America 
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1.  Introduction 

Ultrasound  is  used  extensively  for  differentiation  of 
cysts  from  solid  lesions  in  breast  examinations,  and  it 
is  routinely  used  in  conjunction  with  mammography. 
Ultrasound  can  detect  breast  lesions  a  few  millime¬ 
ters  in  size.1  However,  its  specificity  in  breast  can¬ 
cer  diagnosis  is  not  considered  to  be  high  enough  as  a 
result  of  overlapping  characteristics  of  benign  and 
malignant  lesions.2-3  Optical  imaging  based  on  dif¬ 
fusive  near-infrared  (NIR)  light  has  the  great  poten¬ 
tial  to  differentiate  tumors  from  normal  breast 
tissues  through  determination  of  tissue  parameters, 
such  as  blood  volume,  blood  02  saturation,  tissue 
light  scattering,  water  concentration,  and  the  concen¬ 
tration  and  lifetime  of  exogenous  contrast  agents.4-12 
As  a  potential  diagnostic  tool,  however,  NIR  diffusive 
light  imaging  suffers  from  low  spatial  resolution  and 
lesion  location  uncertainties  because  of  intense  light 
scattering  in  tissue. 
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Most  NIR  imaging  reconstruction  algorithms  are 
based  on  tomographic  inversion  techniques.13-20 
Reconstruction  of  tissue  optical  properties  in  general 
is  underdetermined  and  ill-posed  because  the  total 
number  of  unknown  optical  properties  always  ex¬ 
ceeds  the  number  of  measurements,  and  the  pertur¬ 
bations  produced  by  the  heterogeneities  are  much 
smaller  than  the  background  signals.  In  addition, 
the  inversion  reconstruction  algorithms  are  sensitive 
to  measurement  noise  and  model  errors. 

Our  group  and  others  have  introduced  a  novel  hy¬ 
brid  imaging  method  that  combines  the  complemen¬ 
tary  features  of  ultrasound  and  NIR  diffusive  light 
imaging.21^25  The  hybrid  imaging  obtains  coregis¬ 
tered  ultrasound  and  NIR  diffusive  light  images 
through  simultaneous  deployment  of  an  ultrasound 
array  and  NIR  source- detector  fibers  on  the  same 
probe.21*22*24  Coregistration  permits  joint  evalua¬ 
tion  of  acoustic  and  optical  properties  of  breast  le¬ 
sions  and  enables  use  of  lesion  morphology  provided 
by  high-resolution  ultrasound  to  improve  the  lesion 
optical  property  estimate.  With  the  a  priori  knowl¬ 
edge  of  lesion  location  and  shape  provided  by  coreg¬ 
istered  ultrasound,  NIR  imaging  reconstruction  can 
be  localized  within  specified  three-dimensional  (3-D) 
regions.  As  a  result,  the  reconstruction  is  over  deter¬ 
mined  because  the  total  number  of  unknown  optical 
properties  is  reduced  significantly.  In  addition,  the 
reconstruction  is  less  sensitive  to  noise  because  the 
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convergence  can  be  achieved  with  a  small  number  of 
iterations. 

The  clinical  use  of  the  combined  diagnosis  relies  on 
the  coregistration  of  both  ultrasound  and  NIR  sen¬ 
sors  at  the  probe  level  Conventional  ultrasound 
pulse-echo  imaging  requires  that  an  imaging  trans¬ 
ducer  be  located  on  top  of  the  target,  whereas  NIR 
diffusive  light  imaging  is  feasible  when  the  optical 
source  and  detector  fibers  are  distributed  at  the  pe¬ 
riphery  of  the  ultrasound  transducer.  However,  the 
effects  of  missing  optical  sources  in  the  middle  of*  the 
combined  probe  on  the  accuracy  of  the  reconstructed 
optical  properties  have  to  be  evaluated.  In  addition, 
the  improvements  of  reconstructed  optical  properties 
with  the  guidance  of  the  coregistered  ultrasound 
need  to  be  quantitatively  assessed.  Furthermore, 
real-time  data  acquisition  is  necessary  to  avoid  errors 
in  coregistration  caused  by  patient  motion  during  the 
clinical  experiments.  In  this  paper  we  report  our 
experimental  results  on  the  optimal  probe  configura¬ 
tion,  and  we  quantify  the  improvements  on  recon¬ 
structed  optical  properties  using  a  combined  probe. 
We  also  demonstrate  simultaneous  combined  imag¬ 
ing  with  a  near-real-time  imager. 


age,  and  the  speed  of  reconstruction  has  been 
improved  by  an  order  of  magnitude.  In  this  paper 
we  furthermore  demonstrate  that,  with  the  3-D  tar¬ 
get  distribution  provided  by  coregistered  ultrasound, 
the  accuracy  of  reconstructed  jxa  and  the  reconstruc¬ 
tion  speed  can  be  further  improved. 

To  solve  the  unknown  optical  properties  of  Eq.  (1), 
we  used  the  total  least-squares  (TLS)  method26*27  to 
iteratively  invert  Eq.  (1).  The  TLS  method  performs 
better  than  other  least-squares  when  the  measure¬ 
ment  data  are  subject  to  noise  and  the  linear  operator 
W  contains  errors.  We  found  that  the  TLS  method 
provides  more  accurate  reconstructed  optical  proper¬ 
ties  than  other  least-squares  methods,  and  we 
adopted  TLS  in  solving  inverse  problems.  It  has 
been  shown  in  Ref.  28  that  the  TLS  minimization  is 
equivalent  to  the  following  minimization  problem: 


min 


ll^sd  -  wxf 
IWI2 + 1 


where  X  represents  unknown  optical  properties. 
The  conjugate  gradient  technique  was  employed  to 
iteratively  solve  Eq.  (2). 


2.  Near-Infrared  Diffusive  Wave  Imaging 

We  used  the  Bom  approximation  to  relate  the  scat¬ 
tered  field  UJ  (r,  a))  measured  at  the  probe  surface  to 
absorption  variations  in  each  volume  element  within 
the  sample.  In  the  Bom  approximation,  the  scat¬ 
tered  wave  originated  from  a  source  at  rsi,  and  mea¬ 
sured  at  rdi  it  can  be  related  to  the  medium 
absorption  heterogeneity  A|xa(r^)  at  rvJ  by 

[ffsdWl  =  (1) 

where  M  is  the  total  number  of  source- detector  pairs, 
N  is  the  total  number  of  imaging  voxels,  and  WLj  = 
G(rvJ,  rdi)  <a)Uinc(rvj,  rsi,  <o)vAr V3/D  is  the  weight  ma¬ 
trix  given  in  Ref.  19.  G(rvJ,  rdi,  co)  and  Uiac( rvJ,  rsi,  co) 
are  a  Green’s  function  and  incident  wave,  respec¬ 
tively.  co  is  the  modulation  frequency  and  D  is  the 
average  or  background  diffusion  coefficient,  which  is 
the  average  value  over  the  background  or  whole  tis¬ 
sue. 

With  M  measurements  obtained  from  all  possible 
source- detector  pairs  in  the  planar  array,  we  can 
solve  N  unknowns  of  jia  by  inverting  the  above  matrix 
equation.  In  general,  the  perturbation  Eq.  (1)  is  un¬ 
derdetermined  {M  <  N)  and  ill-posed. 

NIR  imaging  by  itself  generally  has  poor  depth 
discrimination.  However,  ultrasound  is  excellent  in 
providing  accurate  target  depth.  Once  the  target 
depth  is  available  from  coregistered  ultrasound,  we 
can  set  Aga  of  a  nontarget  depth  equal  to  zero.  This 
implies  that  all  the  measured  perturbations  originate 
from  the  particular  depth  that  contains  the  target. 
Because  the  number  of  unknowns  is  reduced  signifi¬ 
cantly,  the  reconstruction  converges  very  fast.  In 
Ref.  23  we  reported  that,  with  a  priori  target  depth 
provided  by  ultrasound,  the  accuracy  of  the  recon¬ 
structed  pa  has  been  improved  by  15-30%  on  aver- 


3.  Methods 

A.  Combined  Probe  and  Imaging  Geometry 

There  are  four  basic  requirements  to  guide  the  design 
of  the  combined  probe.  First,  reflection  geometry  is 
preferred  because  a  conventional  ultrasound  scan  is 
performed  with  this  geometry.  Second,  an  ultra¬ 
sound  array  needs  to  occupy  the  center  of  the  com¬ 
bined  probe  for  coherent  imaging.  Third,  NIR 
sources  and  detectors  have  to  be  distributed  at  the 
periphery.  Because  photon  propagation  distribu¬ 
tion  exhibits  a  banana  shape,  imaging  of  the  tissue 
volume  underneath  the  probe  is  feasible  even 
through  there  are  no  sources  and  detectors  deployed 
in  the  central  portion  of  the  probe.  Fourth,  the  min¬ 
imum  source- detector  separation  should  be  larger 
than  1  cm  for  the  diffusion  approximation  to  be  valid, 
and  the  maximum  separation  should  be  ~8-9  cm  to 
effectively  probe  depths  of  3-4  cm. 

On  the  basis  of  these  requirements  we  deployed  12 
dual-wavelength  optical  source  fibers  and  8  detector 
fibers  over  a  9  cm  X  9  cm  probe  area  (see  Fig.  1). 
The  minimum  and  maximum  source- detector  sepa¬ 
rations  in  the  configuration  are  1.4  and  8  cm,  respec¬ 
tively.  To  study  the  effect  of  the  central  optical  hole 
on  the  accuracy  of  the  reconstructed  optical  proper¬ 
ties,  we  compared  the  reconstruction  results  with  an 
extra  center  source  and  without  the  center  source. 
The  configuration  without  the  center  source  corre¬ 
sponds  to  a  2  cm  X  2  cm  hole  area.  We  further 
moved  the  noncenter  12  sources  and  8  detectors  to¬ 
ward  periphery  by  leaving  a  3  cm  X  3  cm  hole  area  in 
the  middle.  Figure  2  shows  the  picture  of  a  com¬ 
bined  probe  with  the  3  cm  X  3  cm  central  area  occu¬ 
pied  by  an  ultrasound  array.  The  ultrasound  array 
consists  of  64  elements  made  of  1.5-mm-diameter  pi¬ 
ezoelectric  transducers  (Valpey  Fisher  Inc).  The 
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Fig.  1.  Schematic  arrangement  of  NIR  source  and  detector  fibers 
on  the  probe.  Small  solid  circles  are  the  source  fibers  and  larger 
solid  cycles  are  the  detector  fibers. 


transducers  are  deployed  in  a  rectangular  matrix 
with  4-mm  spacing  in  both  x  and  y  directions.  The 
center  frequency  of  the  transducer  is  6  MHz  and  the 
bandwidth  is  40%.  The  transducers  are  made  from 
the  same  piece  of  piezoelectric  transducer  material. 
Therefore  the  gain  difference  among  different  trans¬ 
ducers  is  less  than  3  dB.  The  12  dual-wavelength 
optical  laser  diode  sources  (760  and  830  nm)  and  8 
photomultiplier  tube  (PMT)  detectors  are  coupled  to 
the  probe  through  optical  fibers,  which  are  deployed 
at  the  periphery  of  the  two-dimensional  (2-D)  ultra¬ 
sound  array.  This  hybrid  array  deployment  compro¬ 
mises  ultrasound  coherent  imaging  and  NIR  diffusive 
light  imaging  characteristics. 

The  9  cm  X  9  cm  x  4  cm  image  volume  underneath 
the  probe  is  discretized  into  voxels  of  size  0.4  cm  X  0.4 
cm  X  1  cm.  There  is  a  trade-off  between  the  accu¬ 
rate  estimation  of  the  weight  matrix  W  and  the  voxel 


Fig.  2.  Picture  of  an  experimental  probe .  An  ultrasound  array  of 
8  X  8  =  64  transducers  occupies  the  central  3  cm  X  3  cm  area,  and 
12  dual-wavelength  source  fibers  and  8  detector  fibers  are  deployed 
at  the  periphery. 


Breast:  patient  in  supine 


Fig.  3.  Schematic  of  the  NIR  frequency- domain  imaging  system. 
The  modulation  frequency  is  140  MHz.  The  12  dual- wavelength 
source  channels  are  switched  on  sequentially  by  a  PC,  and  8  de¬ 
tector  channels  receive  signals  in  parallel.  BPF,  bandpass  filter; 
OSC,  oscillator. 


size.  Because  Wtj  is  a  discrete  approximation  of  the 
integral 

|  G( r„,  rd,  o))[/inc(r„,  r3,  to)  ^  dr?, 

it  is  more  accurate  when  the  voxel  size  is  smaller. 
However,  the  total  number  of  reconstructed  un¬ 
knowns  will  increase  dramatically  with  the  decreas¬ 
ing  voxel  size.  Furthermore,  the  rank  of  W  does  not 
increase  in  the  same  order  as  the  total  number  of 
voxels  when  the  voxel  size  decreases.  This  suggests 
that  neighboring  WtJ s  are  correlated  when  the  voxel 
size  is  smaller,  and  a  further  decrease  in  voxel  size 
will  not  add  more  independent  information  to  the 
weight  matrix.  We  found  that  a  0.4  cm  x  0.4  cm  X 
1  cm  voxel  size  is  a  good  compromise.  Therefore  we 
used  this  voxel  size  in  image  reconstructions  reported 
in  this  paper. 

B.  Experimental  Systems 

i.  Near-Infrared  Imaging  System 
We  constructed  a  NIR  frequency-domain  imaging 
system.  The  block  diagram  of  the  system  is  shown 
in  Fig.  3.  This  system  has  12  dual-wavelength 
source  channels  and  8  parallel  receiving  channels. 
On  the  transmission  part,  12  pairs  of  dual¬ 
wavelength  (780  and  830  nm)  laser  diodes  are  used  as 
light  sources,  and  their  outputs  are  amplitude  mod¬ 
ulated  at  140.000  MHz.  Each  one  of  the  12  optical 
combiners  (OZ  Optics  Inc.)  looks  like  a  Y  adapter, 
guiding  the  emission  of  two  diodes  of  different  wave- 
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lengths  through  the  same  thin  optical  fiber  (approx¬ 
imately  0.2  mm  in  diameter).  To  reduce  noise  and 
interference,  an  individual  driving  circuit  is  built  for 
each  diode.  As  a  laser  diode  works  in  series,  a  con¬ 
trol  board  that  interprets  instructions  from  a  PC  is 
used  to  coordinate  operations  of  associated  compo¬ 
nents.  When  a  single  transmission  channel  is  se¬ 
lected,  it  turns  on  the  corresponding  driving  circuit  so 
that  a  dc  driving  current  can  be  set  up  for  the  diode. 
At  the  same  time,  a  selected  signal  is  sent  to  a  rf 
switching  unit,  which  distributes  a  rf  signal  to  the 
right  channel  to  modulate  the  optical  output.  On 
the  reception  part,  eight  PMTs  are  employed  to  detect 
diffusely  reflected  light  from  turbid  media.  Each 
PMT  is  housed  in  a  sealed  aluminum  box,  shielding 
both  environment  lights  and  electromagnetic  fields, 
and  an  optical  fiber  (3  mm  in  diameter)  couples  NIR 
light  from  the  detection  point  to  the  reception  window 
of  the  PMT.  The  electrical  signal  converted  from  the 
optical  input  is  generally  weak  and  rather  high  in 
frequency,  so  high-gain  amplification  and  frequency 
transform  are  necessary  before  it  can  be  sampled  by 
an  analog-to-digital  (A/D)  board  inside  the  PC.  We 
built  eight  parallel  heterodyne  amplification  chan¬ 
nels  to  measure  the  response  of  all  detectors  simul¬ 
taneously,  which  reduces  the  data-acquisition  time. 
Each  amplification  channel  consists  of  a  rf  amplifier 
(40  dB),  a  mixer  in  which  the  rf  signal  (OSC1, 140.000 
MHz)  is  mixed  with  a  local  oscillator  (OSC2, 140.020 
MHz),  a  bandpass  filter  centered  at  20  kHz,  and  a 
low-frequency  amplifier  of  30  dB.  The  heterodyned 
two-stage  amplification  scheme  helps  suppress  wide¬ 
band  noises  efficiently.  We  also  generated  a  refer¬ 
ence  signal  of  20  kHz  by  directly  mixing  OSCl  and 
OSC2,  which  is  necessary  for  retrieving  phase  shifts. 
Eight  detection  signals  and  one  reference  are  sam¬ 
pled,  converted,  and  acquired  into  the  PC  simulta¬ 
neously,  in  which  the  Hilbert  transform  is  used  to 
compute  the  amplitude  and  phase  of  each  channel. 
The  entire  data  acquisition  takes  less  than  1  min, 
which  is  fast  enough  to  acquire  data  from  patients. 

One  of  the  challenges  encountered  in  the  design  of 
a  NIR  imaging  system  is  the  huge  dynamic  range  of 
signals  received  at  various  source- detector  dis¬ 
tances.  For  example,  for  a  semi-infinite  phantom 
made  of  0.5%  Intralipid  solution,  the  amplitude  mea¬ 
sured  at  1  cm  away  from  a  source  is  approximately 
5000  times  larger  than  that  at  8-cm  separation.  In 
addition,  the  perturbation  that  is  due  to  an  embedded 
heterogeneity  with  optical  properties  similar  to  a  tu¬ 
mor  is  normally  a  few  percent  of  the  background 
signal.  As  a  result,  a  reflection-mode  NIR  imaging 
system  should  have  at  least  a  120-dB  dynamic  range 
to  probe  a  target  up  to  4  cm  in  depth.  It  is  hard  to 
build  amplifiers  that  work  linearly  over  such  a  wide 
dynamic  range.  We  overcome  this  difficulty  by  im¬ 
plementing  two-level  source  outputs.  The  dc  output 
of  a  laser  diode  is  controlled  when  its  feedback  loop  is 
adjusted,  whereas  the  rf  signal  is  switched  simulta¬ 
neously  by  a  two-step  attenuator  (no  attenuation  or 
30-dB  attenuation).  When  the  source  and  detector 
are  close  to  each  other,  the  source  is  controlled  to 


yield  a  low-level  output.  When  the  separation  be¬ 
comes  larger,  a  30-dB  higher  output  level  should  be 
used.  With  this  two-level  source  scheme,  our  system 
achieved  fairly  good  linearity  over  a  wide  range  of 
source- detector  separations  (from  1.5  to  8  cm). 

Because  the  parameters  of  an  individual  laser  di¬ 
ode  or  a  PMT  vary  considerably  from  one  to  another, 
we  have  to  calibrate  the  gain  and  phase  shift  for  each 
channel.  A  set  of  measurements  obtained  from  all 
source- detector  pairs  placed  on  the  boundary  of  a 
homogeneous  medium  is 

Aup,  a  =  1,  2,  ...  ,  77i,  (J  —  1,  2,  . . .  ,  /i. 

Here,  amplitude  AaP  and  phase  $aP  are  related  to 
source  a  and  detector  [J,  and  m  and  n  are  the  total 
number  of  sources  and  detectors,  respectively. 
From  the  diffusion  theory,  we  can  obtain  the  follow¬ 
ing  set  of  equations7: 


Aa&  =  Js(a)Jd(p) 


exp(-/e1pa3) 

2 

Pa{3 


4>a[J  =  <Ps(<*)  +  9d(P)  +  *rPa&, 

in  which  Is{ a)  and  <pfi(a)  are  the  relative  gain  and 
phase  delay  associated  with  source  channel  a,  Jd((3) 
and  <pd(p)  are  similar  quantities  associated  with  de¬ 
tector  channel  p,  paP  is  the  corresponding  separation, 
and  kr  +  jkt  is  the  complex  wave  number.  We  obtain 
the  following  set  of  linear  equations  by  taking  a  log¬ 
arithm  of  the  above  equations  related  to  amplitude: 


log(pap2Aa|3)  =  log(/s(a)]  +  logf/dO)]  -  &,Pap, 

4>afJ  =  9.(0)  +  <Pd(P)  +  ArPcp.  (3) 

Although  the  optical  properties  of  the  calibration  me¬ 
dium  are  known  in  advance,  we  leave  the  wave  num¬ 
ber  as  a  variable  and  use  fitted  kr  and  ki  to  calculate 
the  background  scattering  and  absorption  coeffi¬ 
cients.  We  verified  our  calibration  method  by  com¬ 
paring  the  best  fitted  with  real  values.  The 
results  of  our  using  0.5- 0.8%  Intralipid  solutions  al¬ 
ways  yielded  scattering  and  absorption  coefficients 
with  a  rather  good  accuracy.  With  the  two  unknown 
wave  numbers  included,  the  total  number  of  un¬ 
knowns  is  2(m  -h7i  +  2),  which  is  generally  far 
smaller  than  the  number  of  measurements  m  X  n. 
Consequently,  Eq.  (3)  is  overdetermined.  We  can 
solve  all  J»,  Jd(p),  <p»,  and  cp^(p)  terms  as  well  as 
two  unknown  wave  numbers  in  a  least-squares  sense. 
Then  all  measurements  can  be  calibrated  accord¬ 
ingly.  The  results  of  amplitude  AaP  =  exp(-&tpaP)/ 
paP2  and  phase  4>ap  =  fcrpaP  after  calibration  are 
shown  in  Fig.  4.  As  one  can  see,  the  calibrated  am¬ 
plitude  (log  paP2Aa(J)  and  the  phase  from  various 
source- detector  pairs  change  linearly  with  distance. 


2.  Ultrasound  System 

The  ultrasound  system  diagram  is  shown  in  Fig.  5, 
and  the  system  consists  of  64  parallel  transmission 
and  receiving  channels.  Each  transmission  circuit 
can  generate  a  high-voltage  pulse  of  200-ns  duration 
(6  MHz)  with  125  V  peak  to  peak  to  the  connected 
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Fig.  4.  (a)  Log(pQ32AQp)  versus  distancev  pa(3  after  calibration, 
(b)  Phase  c}>fx(j  versus  distance  paP  after  calibration. 


transducer.  Each  receiving  circuit  has  two-stage 
amplifiers  followed  by  an  A/D  converter  with  40-MHz 
sampling  frequency.  The  amplifier  gain  can  be  con¬ 
trolled  based  on  the  target  strength.  A  group  of 
transmission  channels  can  be  addressed  simulta¬ 
neously  to  transmit  pulses  from  neighbor  transduc¬ 
ers  with  specified  delays  and  therefore  to  focus  the 
transmission  beam.  The  re  tuned  signals  can  be  re¬ 
ceived  simultaneously  by  a  group  of  transducers,  and 
the  signals  can  be  summed  with  specified  delays  to 
form  a  receiving  beam. 

The  data-acquisition  procedure  is  the  following. 
The  first  9-element  neighbor  subarray  (dashed  rect¬ 
angle  in  Fig.  6)  from  the  64-element  transducer  array 
and  the  corresponding  channels  are  chosen,  and  then 
the  transmission  delay  profiles  are  generated  in  the 
computer  according  to  the  prespecified  focal  depth. 
The  delay  profile  data  are  transferred  to  the  64- 
channel  delay  profile  generator,  which  triggers  the  64 
high-voltage  pulsers  as  well  as  the  receiving  chan¬ 
nels.  The  returned  ultrasound  signals  are  amplified 


Fig.  5.  Schematic  of  our  ultrasound  scanner.  We  connected  64 
ultrasound  transducers  to  64  parallel  transmission  and  reception 
channels.  The  transmission  part  consists  of  64  high-voltage 
pulsers,  which  can  be  controlled  by  computer-generated  delay  pro¬ 
files.  The  reception  part  consists  of  64  two-stage  amplifiers  and 
A/D  converters.  CH,  channel. 


by  two-stage  amplifiers  and  sampled  by  A/D  convert¬ 
ers.  The  data  are  buffered  in  the  memories  and  are 
read  by  the  computer  after  the  entire  data  acquisition 
is  completed.  The  second  subarray  (solid  rectangle 
in  Fig.  6)  is  chosen  and  the  same  data-acquisition 
process  is  repeated.  A  total  of  64  subarrays  is  used 
in  the  data  acquisition.  After  the  64-subarray  data 
acquisition  is  completed,  the  data  stored  in  the  mem¬ 
ories  are  read  by  the  computer  for  image  formation. 
The  entire  data  acquisition  and  imaging  display  are 
performed  in  approximately  5  s,  which  is  fast  enough 
for  clinical  experiments.  To  ensure  good  signal-to- 
noise  ratio,  we  perform  all  the  electronics  using 
printed  circuit  boards. 

Figure  7  shows  the  picture  of  the  entire  system  and 
the  combined  probe.  Both  the  NIR  system  (top)  and 
the  ultrasound  system  (bottom)  are  mounted  on  a 
hospital  cart.  The  combined  probe,  which  houses 
the  ultrasound  array  and  the  NIR  source-detector 
fibers,  is  designed  to  be  hand  held  to  scan  patients. 
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© 
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Fig.  6.  Ultrasound  subarray  scanning  configuration. 
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Fig.  7.  Picture  of  our  combined  system.  NIR  system  (top)  and 
ultrasound  system  (bottom)  are  mounted  on  a  hospital  cart. 


C.  Phantoms 

We  used  0. 5-0.6%  Intralipid  solutions  to  mimic  nor¬ 
mal  human  breast  tissues  in  all  experiments,  and  the 
corresponding  reduced  scattering  coefficient  jjls' 
ranges  from  5  to  6  cm"1.  The  Intralipid  is  contained 
in  a  large  fish  tank  to  set  up  approximately  a  semi¬ 
infinite  homogeneous  phantom.  Small  semispheri- 
cal  balls  (1  cm  in  diameter),  made  of  acrylamide  gel,22 
are  inserted  into  Intralipid  to  emulate  lesions  embed¬ 
ded  in  a  breast.  The  reduced  scattering  coefficients 
of  the  gel  phantoms  are  similar  to  that  of  the  back¬ 
ground  medium  (ps'  »  6  cm"1),  and  we  changed  the 
absorption  coefficients  to  different  values  by  adding 
different  concentrations  of  India  ink  to  emulate  high- 
contrast  (jxa  =  0.25-cm"1)  and  low-contrast  (\xa  = 
0,1-cm"1)  lesions.  Ultrasound  scattering  particles 
of  200  fxm  in  diameter  are  added  to  the  gel  phantom 
before  the  gel  is  formed. 

4.  Experimental  Results 

A.  Effects  of  Missing  Optical  Sources  in  the  Combined 
Probe 

A  series  of  experiments  was  conducted  to  estimate 
the  optimal  hole  size.  Three  probe  configurations 
were  investigated:  (a)  no-hole,  (b)  2  cm  X  2  cm  cen¬ 
tral  hole,  and  (c)  3  cm  X  3  cm  hole  probes.  The 
no-hole  probe  was  essentially  the  same  as  case  (b) 
except  that  an  additional  light  source  was  added  in 
the  middle.  Figure  8  shows  reconstructed  NIR  im¬ 
ages  for  on-center  targets  of  high  (|i.a  =  0.25  cm"1,  left 
column)  and  low  contrast  (p.a  =  0.1  cm"1,  right  col¬ 
umn)  located  2.5  cm  deep  inside  the  Intralipid.  The 
fitted  background  jxa  and  |xs  are  0.015  and  5.36  cm"1, 
respectively.  With  the  target  depth  provided  by  ul¬ 
trasound,  we  performed  reconstruction  in  the  target 
layer.  The  centers  of  the  voxels  in  this  layer  were  (x, 
y ,  2.5  cm),  where  x  and  y  were  discrete  spatial  x-y 
coordinates,  and  the  thickness  of  the  layer  was  1  cm. 
For  the  high-contrast  target  case,  there  are  no  im¬ 
portant  differences  in  image  quality  associated  with 
different  probes  [Figs.  8(a)  and  8(c)]  except  that  with 
a  3  cm  X  3  cm  hole.  The  first  row  of  Table  1  provides 


measured  maximum  p.a  values  from  the  correspond¬ 
ing  images.  Because  of  the  low  spatial  resolution  of 
diffusive  imaging,  the  boundaries  of  the  targets  are 
not  well  defined.  The  maximum  value  is  a  better 
estimation  of  reconstructed  target  fxa.  From  no  hole 
to  2  cm  x  2  cm,  the  reconstructed  maximum  jji0  de¬ 
creases  slowly.  But  for  3  cm  X  3  cm,  the  maximum 
\ia  drops  suddenly  to  0.104  cm"1,  which  is  less  than 
half  of  the  original  value.  Another  imaging  param¬ 
eter  we  measured  is  the  full  width  at  half-maximum 
(FWHM)  of  the  corresponding  images.  Because  the 
image  lobes  were  elliptical  in  general,  we  measured 
the  widths  of  longer  and  shorter  axes  and  used  the 
geometric  mean  to  estimate  the  FWHM.  The  re¬ 
sults  are  shown  in  Table  1,  and  the  FWHM  almost 
increases  with  the  hole  size.  We  also  measured  the 
image  artifact  level,  which  was  defined  as  the  ratio  of 
the  peak  artifact  to  the  maximum  strength  of  the 
image  lobe  and  is  given  in  decibels.  The  results  are 
shown  in  Table  1.  No  artifacts  were  observed  in  the 
images  of  no-hole  and  2  cm  X  2  cm  hole  probes. 
However,  the  peak  artifact  level  at  the  -  14.3-dB  level 
was  measured  in  the  image  of  the  3  cm  X  3  cm  hole 
probe.  When  the  contrast  was  low,  the  recon¬ 
structed  maximum  absorption  coefficients  and  mea¬ 
sured  FWHMs  were  essentially  the  same  for  the 
no-hole  and  2  cm  X  2  cm  hole  probes.  However,  the 
reconstructed  maximum  value  dropped  to  60%  of  the 
true  value  for  the  3  cm  X  3  cm  probe.  The  artifact 
levels  measured  in  the  images  of  three  probe  config¬ 
urations  were  similar  and  were  worse  than  the  high- 
contrast  case.  The  image  artifacts  are  related  to  the 
reconstruction  algorithm.  When  the  target  contrast 
is  weak  or  the  signal-to-noise  ratio  is  low,  the  inver¬ 
sion  algorithm  produces  artifacts  around  the  edges  of 
the  images. 

For  shallow  targets  (here  we  set  the  target  depth  to 
be  1.5  cm)  the  NIR  system  has  a  relatively  poorer 
performance.  This  is  due  to  less  source- detector 
pairs  experiencing  the  existence  of  a  shallow  absorber. 
As  shown  in  Fig.  9,  image  artifacts  are  obviously  worse 
compared  with  Fig.  8.  However,  the  conclusion  about 
the  hole  size  of  the  probe  remains  true.  Table  2  lists 
all  the  measured  imaging  parameters  obtained  from 
three  probe  configurations.  Although  a  3  cm  X  3  cm 
hole  is  somewhat  too  big  to  obtain  good  enough  results, 
the  optimal  hole  size  is  in  the  neighborhood  of  2  cm  X 
2  cm.  This  optimal  size  is  approximately  the  size  of 
commercial  ultrasound  transducers. 

In  the  above  studies,  we  used  the  iteration  number 
obtained  from  the  no-hole  configuration  for  the  rest  of 
the  configurations.  Ideally,  the  iteration  should 
stop  when  the  object  function  [see  Eq.  (2)]  or  the  error 
performance  surface  reaches  the  noise  floor.  How¬ 
ever,  system  noise,  particularly  coherent  noise,  was 
difficult  to  estimate  from  experimental  data.  In  gen¬ 
eral,  we  found  that  the  reconstructed  values  were 
closer  to  true  values  when  the  object  function  reached 
approximately  5-15%  of  the  initial  value  (total  en¬ 
ergy  in  the  measurements).  Therefore  we  used  this 
criterion  (—10%  of  the  initial  value)  for  the  no-hole 
configuration.  Because  the  signal-to-noise  ratio  of 
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Fig.  8.  Reconstructed  NIR  images  of  deeper  targets^ (2.5  cm  in  depth,  1  cm  in  diameter,  and  the  fitted  background  p,0  and  \ia  are  0.015 
and  5.36  cm”1,  respectively).  The  left  column  corresponds  to  images  of  a  high-contrast  target  (\x.n  =  0.25  cm”1)  obtained  from  different 
probe  configurations,  and  the  right  column  corresponds  to  images  of  a  low-contrast  target  (p,a  =  0.1  cm*"1).  Each  row  is  related  to  a  specific 
hole  size:  (a)  and  (b)  no  hole,  (c)  and  (d)  2  cm  X  JSj'cm,  (e)  and  (f)  3  cm  x  3  cm. 


the  data  decreased  with  the  increase  in  hole  size,  we 
could  not  find  consistent  criterion  for  both  no-hole 
and  hole  data.  Therefore  we  used  the  same  iteration 
number  obtained  from  the  no-hole  case  for  the,  hole 
configurations,  and  the  comparison  was  based  on  the 
same  iteration  number.  .  ?rf:  * 

B.  Ultrasound-Guided  Near-Infrared  Imaging, 

Three-dimensional  ultrasound  images  can’, provide 
3-D  distributions  of  targets.  With  the  a  priorhtarget 


depth  information,  the  optical  reconstruction  can  be 
improved  significantly.  An  example  is  given  in  Fig. 
10.  The  target  again  was  a  1-cm-diameter  gel  ball  of 
low  (p.a  =  0.1-cnT1)  optical  contrast  and  was  embed¬ 
ded  at  approximately  (0,  0,  2.5  cm)  inside  the  In- 
tralipid  medium.  The  fitted  background  and  \is' 
are  0.02  and  5.08  cm”1,  respectively.  The  combined 
iprobe  shown  in  Fig.  2  was  used  to  obtain  the  ultra¬ 
sound  and  NIR  data  simultaneously.  Figure  10(a) 
shows  an  A-scan  line  of  a  returned  ultrasound  echo 
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Table  1.  Parameters  of  Reconstructed  Images  for  Deep  High-Contrast 
=  0.25-cm"1)  and  Low-Contrast  (pa  =  0.1 -cm-1)  Targets* 


* 

Probe  Type 

Parameter 

No  Hole 

2  cm  X  2  cm 

3  cm  X  3  cm 

High  contrast 

M'a(max) 

(cm"1) 

0.251 

0.234 

0.104 

FWHMC 

(cm) 

1.85 

1.91 

2.44 

Artifacts 

(dB) 

Low  contrast 

Background 

(-22) 

Background 

-14.3 

M'a(mox) 

(cm"1) 

0.105 

0.111 

.0.064 

FWHM 

(cm) 

2.02 

1.83 

2.16 

Artifacts 

(dB) 

-6,90 

-8.10 

-5.65 

“The  fitted  background  p.a  and  p./  are  0.015  and  5.36  cm"1, 
respectively. 

6jia(max)  is  the  measured  maximum  value  of  the  reconstructed 
absorption  coefficient  map. 

CFWHM  is  defined  as  the  geometric  mean  of  the  widths  mea¬ 
sured  at  longer  and  shorter  axes  of  the  elliptical  image  lobe. 


signal  received  by  one  ultrasound  transducer  located 
on  top  of  the  target.  As  acoustic  scatters  were  uni¬ 
formly  distributed  in  the  target,  signals  were  re¬ 
flected  from  inside  the  target  as  well  as  from  the 
surfaces.  The  reflected  signals  from  the  front  and 
back  surfaces  of  the  gel  ball  can  be  clearly  identified 
in  the  echo  signal.  On  the  basis  of  the  target  depth, 
we  reconstructed  the  optical  absorption  coefficient  at 
the  target  depth  only  (1  cm  in  thickness)  by  setting 
the  purturbations  from  the  other  depths  equal  to 
zero.  We  also  performed  3-D  optical-only  recon¬ 
struction.  Figure  10(b)  shows  the  reconstructed  ab¬ 
sorption  image  from  a  3-D  optical-only  reconstruction 
[layer  three  of  voxel  coordinates  (x,  y,  2.5  cm)  and  1 
cm  thick],  whereas  Fig.  10(c)  shows  the  reconstructed 
image  of  the  corresponding  target  from  ultrasound- 
guided  reconstruction.  For  optical-only  reconstruc¬ 
tion,  the  algorithm  did  not  converge  to  a  localized 
spatial  region,  and  the  image  contrast  was  poor. 
The  measured  maximum  absorption  coefficient  was 
0.088  cm*” \  which  was  close  to  the  true  value.  How¬ 
ever,  the  measured  spatial  location  of  the  maximum 
value  was  (— 1.6,  - 1.2  cm),  which  was  too  far  from  the 
true  target  location.  With  the  a  priori  target  depth, 
the  reconstruction  performed  at  the  target  layer  can 
localize  the  target  to  the  correct  spatial  position. 
The  measured  maximum  absorption  coefficient  was 
0.12  cm-1  and  its  location  was  (0,  0.4  cm),  which  was 
very  close  to  the  true  target  location.  This  example 
demonstrates  that  a  priori  target  depth  can  signifi¬ 
cantly  improve  the  reconstruction  accuracy  and  tar¬ 
get  localization. 

In  addition  to  use  of  a  priori  target  depth  informa¬ 
tion,  we  can  also  use  the  target  spatial  distribution 
provided  by  ultrasound  to  guide  the  reconstruction. 
We  performed  a  set  of  experiments  with  two  targets 


located  at  2.5  cm  in  depth  inside  the  Intralipid. 
Each  target  is  a  1-cm3  gel  cube  containing  ultrasound 
scatters.  For  optical  properties,  they  both  could  be 
high  contrast  (pia  ==  0.25  cm"1)  or  low  contrast  (juia  = 
0,1  cm"1),  but  had  the  same  reduced  scattering  coef¬ 
ficient  as  the  background.  The  fitted  background  \x.a 
and  jxs'  are  0.017  and  4.90  cm"1,  respectively.  One 
target  was  centered  approximately  at  (-1.0,  —1.0, 2.5 
cm),  whereas  the  other  was  at  (1.0, 1.0, 2.5  cm).  The 
distance  between  the  centers  of  the  two  targets  was 
2.8  cm. 

Figure  11(a)  is  the  ultrasound  image  of  two  high- 
contrast  targets.  As  the  field  of  view  of  the  ultra¬ 
sound  system  was  nearly  a  3  cm  X  3  cm  square,  these 
two  targets  appeared  at  diagonal  comers.  The  mea¬ 
sured  peak  positions  of  the  two  targets  were  (“0.6, 
“1.0  cm)  and  (1.0,  1.0  cm),  which  differed  from  the 
true  target  locations  by  only  one  voxel.  The  low  con¬ 
trast  of  the  ultrasound  image  is  related  to  the  speckle 
noise.  Because  our  ultrasound  array  is  sparse,  the 
imaging  quality  is  not  state  of  the  art  (see  more  dis¬ 
cussion  in  Section  5).  The  NIR  image  of  these  tar¬ 
gets  was  obtained  simultaneously  and  is  shown  in 
Figure  11(b).  We  performed  the  reconstruction  at 
the  target  layer  by  taking  advantage  of  target  depth 
information.  A  total  of  123  iterations  was  used  to 
obtain  Fig.  11(b).  The  measured  peak  positions  of 
the  two  targets  were  (-1.4,  -1.0  cm)  and  (0.6,  0.6 
cm),  which  were  one  voxel  off  from  the  true  target 
locations  (“1.0,  -1.0  cm)  and  (1.0,  1.0  cm),  respec¬ 
tively.  The  corresponding  reconstructed  absorption 
coefficients  were  0.242  and  0.251  cm*”1,  which  were 
close  to  the  true  values.  However,  the  two  targets 
were  almost  connected  to  each  other,  and  their  spa¬ 
tial  localization  was  poor.  For  low-contrast  targets, 
the  ultrasound  image  is  shown  in  Fig.  11(c),  and  the 
measured  peak  locations  of  the  two  targets  were 
(-1.0,  -0.6  cm)  and  (0.6, 1.0  cm),  which  differed  from 
the  true  target  locations  by  only  one  voxel.  The  cor¬ 
responding  NIR  image  is  shown  in  Fig.  11(d),  and  the 
measured  peak  locations  of  the  two  targets  were 
(-2.2,  -1.0  cm)  and  (0.6,  1.0  cm).  The  left  target 
was  off  the  true  location  by  three  voxels.  The  corre¬ 
sponding  reconstructed  absorption  coefficients  were 
0.063  and  0.1004  cm"1  at  87  iteration  steps.  As  one 
can  see,  the  target  shape  and  localization  were  poorer 
than  those  in  the  high-contrast  case.  In  addition,  an 
artifact  appeared  at  the  edge  of  the  image. 

From  the  coregistered  ultrasound  images,  we  ob¬ 
tained  spatial  distributions  of  the  two  targets  and 
specified  target  regions.  Figures  12(a)  and  12(c) 
show  the  -6-dB  contour  plots  of  Figs.  11(a)  and  11(c). 
Applying  the  same  reconstruction  scheme  to  these 
specific  regions,  we  obtained  Figs.  12(b)  and  12(c)  in 
one  iteration.  The  reconstructed  absorption  coeffi¬ 
cients  were  0.2357  and  0.219  cm"1  for  the  two  high- 
contrast  target  cases  and  0.123  and  0.131  cm"1  for 
the  low-contrast  case.  We  can  see  much  better  im¬ 
provement  in  the  low-contrast  target  case  when  we 
compare  Fig.  12(d)  with  Fig.  11(d).  This  example 
demonstrates  that,  when  the  targets  are  visible  in 
ultrasound  images,  their  morphology  information 
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(e)  (0 

Fig.  9.  Reconstructed  NIR  images  for  shallow  targets  (1.5  cm  in  depth,  1  cm  in  diameter,  and  the  fitted  background  and  \x.,'  are  0.015 
and  5.36  cm-1,  respectively).  The  left  column  corresponds  to  images  of  a  high-contrast  target  (pa  =  0.25  cm  l),  and  the  right  column 
corresponds  to  images  of  a  low-contrast  target  (p.0  =  0.1  cm'1).  Each  row  is  related  to  a  specific  hole  size:  (a)  and  (b)  no  hole,  (c)  and 
(d)  2  cm  x  2  cm,  (e)  and  (f)  3  cm  x  3  cm. 


provided  by  ultrasound  can  be  used  to  guide  the  op¬ 
tical  reconstruction  in  the  specified  regions. 

The  result  regarding  the  iteration  step  is  signifi¬ 
cant.  As  we  discussed  above,  there  is  no  known 
stopping  criterion  to  terminate  the  iteration  because 
it  is  difficult  to  estimate  the  noise  level  in  the  mea¬ 
surements.  With  the  a  priori  target  depth  and  spa¬ 
tial  distribution  provided  by  coregistered  ultrasound, 
we  can  obtain  an  accurate  optical  absorption  coeffi¬ 
cient  in  one  iteration.  Therefore  no  stopping  crite¬ 


rion  is  needed  for  the  inversion  algorithms. 
However,  this  result  will  need  to  be  further  evaluated 
with  more  samples  of  different  contrasts. 

5.  Discussion 

Commercial  ultrasound  scanners  use  one¬ 
dimensional  probes  that  provide  2-D  images  of  x-z 
views  of  the  targets,  where  x  and  z  are  the  spatial  and 
propagation  dimensions,  respectively.  Such  x-z  im¬ 
ages  cannot  coregister  with  NIR  images,  which  are 
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Table  2.  Parameters  of  Reconstructed  Images  for  Shallow 
•  High-Contrast  (|ia  =  0.25-cm-1)  and  Low-Contrast  (|ia  =  0.1-cm-1) 
Targets” 


Probe  Type 

Parameter 

No  Hole 

2  cm  X  2  cm 

3  cm  x  3  cm 

High  contrast 

^a(iaax) 

0.250 

0.194 

0.118 

(cm-1) 

'  FWHM 

1.32 

1.61 

2.08 

(cm) 

Artifacts 

-7.98 

-12.7 

-9.76 

(dB) 

Low  contrast 

P'o(max) 

0.100 

0.091 

0.042 

(cm”1) 

FWHM 

1.88 

2.11 

3.17 

(cm) 

Artifacts 

-6.25 

-7.44 

-0.65 

(dB) 

“The  fitted  background  and  p.s'  are  0.015  and  5.35  cm-1, 
respectively. 


obtained  from  x-y  views  of  the  targets.  Our  current 
2-D  ultrasound  array  is  capable  of  providing  x—y 
views  of  the  targets,  which  can  be  used  to  coregister 
with  NIR  images.  However,  the  array  is  sparse  and 
therefore  the  image  resolution  is  not  state  of  the  art. 
Nevertheless,  its  spatial  resolution  is  comparable  to 
NIR  imaging  and  can  be  used  to  guide  NIR  image 
reconstruction.  With  3-D  ultrasound  guidance,  only 
one  iteration  is  needed  to  obtain  accurate  absorption 
coefficients.  This  result  is  significant  because  no 
stopping  criterion  is  necessary.  More  studies  with  a 
variety  of  target  contrasts  and  locations  will  be  per¬ 
formed  to  verify  this  result. 

We  purchased  a  2-D  state-of-the-art  ultrasound  ar¬ 
ray  of  1280  transducer  elements  and  we  are  building 
a  multiplexing  unit  for  our  64-channel  electronics. 
In  addition,  the  new  2-D  transducer  size  is  approxi¬ 
mately  2  cm  X  3  cm,  which  is  in  the  neighborhood  of 
the  optimal  hole  size  we  found  through  this  study. 
With  the  new  2-D  ultrasound  transducer,  we  will  be 
able  to  obtain  high-resolution  ultrasound  images  and 
delineate  the  target  boundaries  with  finer  details  for 
optical  reconstruction. 

Ultrasound  contrast  depends  on  lesion  acoustic 
properties,  and  NIR  optical  contrast  is  related  to  le¬ 
sion  optical  properties.  Both  contrasts  exist  in  tu¬ 
mors,  but  the  sensitivities  of  these  two  modalities 
may  be  different.  It  is  possible  that  some  early- 
stage  cancers  have  NIR  contrast  but  are  not  detect¬ 
able  by  ultrasound.  It  would  be  desirable  if  we  could 
obtain  sensitivity  of  optical  imaging  alone.  How¬ 
ever,  light  scattering  is  a  main  problem  that  prevents 
the  accurate  and  reliable  localization  of  lesions.  It  is 
also  possible  that  some  lesions  have  acoustic  contrast 
but  no  NIR  contrast  or  low  NIR  contrast.  Currently, 
ultrasound  is  routinely  used  as  an  adjunct  tool  to 
x-ray  mammography;  the  combined  sensitivity  of 
these  two  modalities  in  breast  cancer  detection  is 


more  than  90%.29  Recently,  ultrasound  has  also 
been  advocated  to  screen  dense  breasts.30  We  antic¬ 
ipate  that  our  combined  imaging  will  add  more  spec¬ 
ificity  to  the  ultrasonically  detected  lesions. 

In  the  reported  phantom  studies,  we  assigned  zero 
perturbations  to  the  regions  where  no  targets  were 
present.  In  clinical  studies,  we  plan  to  segment  the 
ultrasound  images  and  specify  different  tissue  types 
as  well  as  suspicious  regions  in  the  segmented  im¬ 
ages.  We  will  then  reduce  the  reconstructed  optical 
unknowns  by  assigning  unknown  optical  properties 
to  different  tissue  types  as  well  as  to  suspicious  re¬ 
gions.  Finally,  we  will  reconstruct  the  reduced  sets 
of  unknown  optical  properties.  We  expect  a  more 
accurate  estimation  of  reconstructed  optical  proper¬ 
ties  and  fast  convergence  speed,  as  reported  in  this 
paper.  However,  it  is  still  too  early  to  judge  the 
clinical  performance  of  the  combined  method;  further 
clinical  studies  are  needed. 

Probing  regions  of  the  banana-shaped  diffusive 
photons  depend  on  source- detector  separations  and 
measurement  geometry.  For  a  semi-infinite  geome¬ 
try,  the  probing  regions  extend  further  into  the  me¬ 
dium  when  source— detector  separation  increases. 
This  is  why  we  have  multiple  source-detector  pairs  of 
various  separations  to  detect  targets  at  variable 
depths  from  0.5  to  4  cm.  Of  course  it  is  hard  to 
achieve  uniform  sensitivity  in  the  entire  region  of 
interest.  For  example,  a  superficial  target  (~1  cm 
deep)  would  cause  strong  perturbations  when  it  is 
close  to  a  source  or  a  detector,  but  will  result  in  much 
weaker  signals  when  it  is  located  deeper.  Normal¬ 
ization  of  scattering  photon  density  waves  with  re¬ 
spect  to  the  incident  waves  makes  it  possible  for 
reconstruction  algorithms  to  handle  the  huge  dy¬ 
namic  range  of  signals  and  to  detect  a  target  as  deep 
as  4  cm.  This  normalization  procedure  was  applied 
to  the  reconstruction  algorithm  used  to  obtain  the 
reported  images. 

In  this  study,  the  target  absorption  coefficient  was 
reconstructed  from  the  measurements.  Because  the 
target  p,s'  was  similar  to  the  background  p,s',  the 
coupling  between  p,a  and  |xs'  in  our  measurements 
was  negligible.  We  also  performed  experiments 
with  gel  phantom  made  with  Intralipid  with  a  con¬ 
centration  similar  to  that  of  the  background  and  did 
not  observe  perturbation  beyond  the  noise  level. 
Similar  reconstruction  studies  can  be  performed  for 
scattering  coefficients  as  well.  Simultaneous  recon¬ 
struction  of  both  absorption  and  scattering  coeffi¬ 
cients  is  also  possible.  Because  the  eigenvalues  of 
the  absorption  and  scattering  weight  matrices  are 
significantly  different,  good  regulation  schemes  are 
needed  for  simultaneous  reconstruction.  This  sub¬ 
ject  is  one  of  our  topics  for  further  study. 

6.  Summary 

We  have  constructed  a  near-real-time  imager  that 
can  provide  coregistered  ultrasound  and  NIR  images 
simultaneously.  This  new  technique  is  designed  to 
improve  the  specificity  of  breast  cancer  diagnosis. 
Because  the  ultrasound  transducer  needs  to  occupy 
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Fig.  10.  Deep  target  (2.5  cm  in  depth,  1  cm  in  diameter)  of  low  optical  contrast  (p0  =  0.10  cm-1  and  fitted  background  and  p/  are  0.02 
and  5.08  cm~\  respectively),  (a)  A-scan  line  of  the  reflected  ultrasound  pulse-echo  signal  indicating  the  target  depth,  (b)  Absorption 
image  of  the  low-contrast  target  obtained  from  optical-only  reconstruction,  (c)  Ultrasound-guided  reconstruction  at  target  depth. 

the  central  region  of  the  combined  probe,  a  series  of  cm,  essentially  similar  reconstruction  results  as 

experiments  were  conducted  to  investigate  the  effects  those  of  no  missing  optical  sensors  in  the  middle  of 

of  missing  optical  sensors  in  the  middle  of  the  com-  the  combined  probe  can  be  obtained.  This  2  cm  X  2 

bined  probe  on  the  NIR  image  quality.  Our  results  cm  dimension  is  approximately  the  size  of  most  com- 

have  shown  that,  as  long  as  the  central  ultrasound  mercial  ultrasound  phased-array  transducers, 

transducer  area  is  in  the  neighborhood  of  2  cm  x  2  When  the  central  missing  optical  sensor  area  is  in- 
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Fig.  11.  Simultaneously  obtained  ultrasound  and  NIK  abson>tion  ima^  ^hetood  backg-ound  | a,  M d»,  ^  ,;d; 

respectively.  W  Ultrasound  „d  <b) I  NIK  .b»rpb.»  ties,  the  W.  «...  «.ro  located  a.  2.5 

NIR  image  of  two  low-contrast  targets  (target  ^  -  0.10  cm  ).  in  Dorn  mgn  an 

cm  in  depth. 


creased  to  3  cm  X  3  cm,  however,  the  reconstructed 
values  are  obviously  lower  than  real  values.  If  we 
increase  the  iteration  steps,  artifacts  in  the  recon¬ 
structed  images  would  soon  become  dominant. 

With  the  target  3-D  distribution  provided  by  coreg¬ 
istered  ultrasound,  significant  improvements  in  algo¬ 
rithm  convergence  and  reconstruction  speed  were 
achieved.  In  general,  a  priori  target  depth  informa¬ 
tion  guides  the  inversion  algorithm  to  reconstruct  the 
heterogeneities  at  the  correct  spatial  locations  and 
improves  the  reconstruction  speed  by  an  order  of 
magnitude.  In  addition,  the  a  priori  target  spatial 
distribution  can  further  reduce  the  iteration  to  one 
step  and  also  obtain  accurate  optical  absorption  co¬ 
efficients.  Given  the  fact  that  no  known  stopping 
criterion  is  available,  this  result  is  significant  be¬ 
cause  no  iteration  is  needed.  However,  this  result 
will  need  to  be  evaluated  with  more  samples  of  dif¬ 
ferent  contrasts. 

We  acknowledge  the  following  for  their  funding 
support:  the  state  of  Connecticut  (99CT21),  U.S. 
Department  of  Defense  Army  Breast  Cancer  Pro¬ 


gram  (DAMD 17-00-1-0217,  DAMD17-01-1-0216), 

and  Multiple-Dimensional  Technology,  Inc. 


References 

1.  T.  A.  Stavros,  D.  Thickman,  and  C.  Rapp,  “Solid  breast  nod¬ 
ules:  use  of  sonography  to  distinguish  between  benign  and 
malignant  lesions,”  Radiology  196,  123-134  (1995). 

2  G  Rahbar,  A.  C.  Sie,  G.  C.  Hansen,  J.  S.  Prince,  M.  L.  Melany, 

'  H.  Reynolds,  V.  P.  Jackson,  J.  W.  Sayre,  and  L.  W.  Bassett, 
“Benign  versus  malignant  solid  breast  masses:  US  differen¬ 
tiation,”  Radiology  213,  889-894  (1999). 

3.  V.  P.  Jackson,  “The  current  role  of  ultrasonography  m  breast 
imaging,”  Radiol.  Clin.  North  Am.  33, 1161-1170  (1995). 

4  B  Tromberg,  N.  Shah,  R.  Lanning,  A.  Cerussi,  J.  Espinoza,  T. 
Pham,  L.  Svaasand,  and  J.  Butler,  “Non-invasive  m  vivo  char¬ 
acterization  of  breast  tumors  using  photon  migration  spectros¬ 
copy,”  Neoplasia  2,  26-40  (2000). 

5  S.  Fantini,  S.  Walker,  M.  Franceschini,  M.  Kaschke,  P.  Schlag 
and  K.  Moesta,  “Assessment  of  the  size,  position,  and  optical 
properties  of  breast  tumors  in  vivo  by  noninvasive  optical 
methods,”  Appl.  Opt.  37, 1982-1989  (1998). 

6.  S.  M.  Nioka,  M.  Shnall,  M.  Miwa,  S.  Orel,  M.  Haida,  S.  Zhao, 
and  B.  Chance,  “Photon  imaging  of  human  breast  cancer,  Adv. 
Exp.  Med.  Biol.  16, 171-179  (1994). 

7.  R.  M.  Danen,  Y.  Wang,  X.  D.  Li,  W.  S.  Thayer,  and  A.  G.  Yodh, 


6378 


APPLIED  OPTICS  /  Vol.  40,  No.  34  /  1  December  2001 


(C)  (d) 

Fig.  12.  (a)  and  (c)  -6-dB  contour  plots  of  ultrasound  images  shown  in  Figs.  11(a)  and  11(c).  The  outer  contour  is  -6  dB  from  the  peak, 
and  the  contour  spacing  is  1  dB.  (b)  and  (d)  Corresponding  NIR  absorption  maps  reconstructed  in  target  regions  specified  by  ultrasound. 


“Regional  imager  for  low  resolution  functional  imaging  of  the 
brain  with  diffusing  near-infrared  light,”  Photochem.  Photo- 
biol.  67,  33-40  (1998). 

8.  T.  McBride,  B.  W.  Pogue,  E.  Gerety,  S.  Poplack,  U.  Osterberg, 
B.  Pogue,  and  K.  Paulsen,  “Spectroscopic  diffuse  optical  tomog¬ 
raphy  for  the  quantitative  assessment  of  hemoglobin  concen¬ 
tration  and  oxygen  saturation  in  breast  tissue,”  Appl.  Opt.  38, 
5480-5490  (1999). 

9.  M.  A.  Franceschini,  K.  T.  Moesta,  S.  Fantini,  G.  Gaida,  E. 
Gratton,  H.  Jess,  M.  Seeber,  P.  M.  Schlag,  and  M.  Kashke, 
“Frequency-domain  techniques  enhance  optical  mammogra¬ 
phy:  initial  clinical  results,”  Proc.  Natl.  Acad.  Sci.  USA  94, 
6468-6473  (1997). 

10.  J.  B.  Fishkin,  O.  Coquoz,  E.  R.  Anderson,  M.  Brenner,  and  B.  J. 
Tromberg,  “Frequency-domain  photon  migration  measure¬ 
ments  of  normal  and  malignant  tissue  optical  properties  in  a 
human  subject,”  Appl.  Opt.  36,  10-20  (1997). 

11.  T.  L.  Troy,  D.  L.  Page,  and  E.  M.  Sevick-Muraca,  “Optical 
properties  of  normal  and  diseased  breast  tissues:  prognosis 
for  optical  mammography,”  J.  Biomed.  Opt  1, 342-355  (1996). 

12.  R.  J.  Grable,  D.  P.  Rohler,  and  S.  Kla,  “Optical  tomography 
breast  imaging,”  in  Optical  Tomography  and  Spectroscopy  of 
Tissue:  Theory ,  Instrumentation,  Model ,  and  Human  Studies 
//,  B.  Chance  and  R.  Alfano,  eds.,  Proc.  SPIE  2979,  197-210 
(1997). 

13.  Y.  Yao,  Y.  Wang,  Y.  Pei,  W.  Zhu,  and  R.  L.  Barbour, 
“Frequency-domain  optical  imaging  of  absorption  and  scatter¬ 


ing  distributions  by  a  Bom  iterative  method,”  J.  Opt.  Soc.  Am. 
A  14,  325-341  (1997). 

14.  H.  Jiang,  K.  Paulsen,  U.  Osterberg,  B.  Pogue,  and  M.  Patter¬ 
son,  “Optical  image  reconstruction  using  frequency-domain 
data:  simulations  and  experiments,”  J.  Opt.  Soc.  Am.  A.  12, 
253-266  (1995). 

15.  X.  Li,  T.  Durduran,  A.  Yodh,  B.  Chance,  and  D.  N.  Pattanayak, 
“Diffraction  tomography  for  biomedical  imaging  with  diffuse- 
photon  density  waves,”  Opt.  Lett.  22,  573-575  (1997). 

16.  C.  Matson  and  H.  Liu,  “Backpropagation  in  turbid  media,”  J. 
Opt.  Soc.  Am.  A  16,  1254-1265  (1999). 

17.  M.  A.  O'Leary,  “Imaging  with  diffuse  photon  density  waves,” 
•  Ph.D.  dissertation  (University  of  Pennsylvania,  Philadelphia, 

Pa.,  1996). 

18.  K.  Paulsen,  P.  Meaney,  M.  Moskowitz,  and  J.  Sullivan,  Jr.,  “A 
dual  mesh  scheme  for  finite  element  based  reconstruction  al¬ 
gorithms,”  IEEE  Trans.  Med.  Imaging  14,  504-514  (1995). 

19.  S.  Arridge  and  M.  Schweiger,  “Photon-measurement  density 
functions.  Part  I:  Analytical  forms,”  Appl.  Opt.  34,  7395- 
7409  (1995). 

20.  S.  Arridge  and  M.  Schweiger,  “Photon-measurement  density 
functions.  II.  Finite-element-method  calculations,”  Appl. 
Opt.  34,  8026-8037  (1995). 

21.  Q.  Zhu,  T.  Dunrana,  M.  Holboke,  V.  Ntziachristos,  and  A. 
Yodh,  “Imager  that  combines  near-infrared  diffusive  light  and 
ultrasound,”  Opt.  Lett.  24,  1050-1052  (1999). 

22.  Q.  Zhu,  D.  Sullivan,  B.  Chance,  and  T.  Dambro,  “Combined 


1  December  2001  /  Vol.  40,  No.  34  /  APPLIED  OPTICS  6379 


*  ultrasound  and  near  infrared  diffusive  light  imaging,”  IEEE 
Trans.  Ultrason.  Ferroelectr.  Freq.  Control  46,  665-678 
(1999). 

23.  Q.  Zhu,  E.  Conant,  and  B.  Chance,  “Optical  imaging  as  an 
adjunct  to  sonograph  in  differentiating  benign  from  malignant 
breast  lesions,"  J.  Biomed.  Opt.  5(2),  229-236  (2000). 

24.  Q.  Zhu,  N.  G.  Chen,  D.  Q.  Piao,  P.  Y.  Guo,  and  X.  H.  Ding, 
“Design  of  near-infrared  imaging  probe  with  the  assistance  of 
ultrasound  localization,"  Appl.  Opt.  40,  3288-3303  (2001). 

25.  M.  Jholboke,  B.  J.  Tromberg,  X.  Li,  N.  Shah,  J.  Fishkin,  D. 
Kidney,  J.  Butler,  B.  Chance,  and  A.  G.  Yodh,  “Three- 
dimensional  diffuse  optical  mammography  with  ultrasound 
localization  in  human  subject,"  J.  Biomed.  Opt.  5(2),  237-247 
(2000). 

26.  W.  Zhu,  Y.  Wang,  and  J.  Zhang,  “Total  least-squares  recon¬ 


struction  with  wavelets  for  optical  tomography,"  J.  Opt.  Soc. 
Am.  A  15,  2639-2650  (1998). 

27.  P.  C.  Li,  W.  Flax,  E.  S.  Ebbini,  and  M.  O'Donnell,  “Blocked 
element  compensation  in  phased  array  imaging,”  IEEE  Trans. 
Ultrason.  Ferroelectr.  Freq.  Control  40(4),  283-292  (1993). 

28.  G.  H.  Golub,  “Some  modified  matrix  eigenvalue  problems,” 
SIAM  (Soc.  Ind.  Appl.  Math.)  Rev.  15,  318-334  (1973). 

29.  H.  Zonderland,  E.  G.  Coerkamp,  J.  Hermans,  M.  J.  van  de 
Vijver,  and  A.  E.  van  Voorthuisen,  “Diagnosis  of  breast  cancer: 
contribution  of  US  as  an  adjunct  to  mammography,”  Radiology 
213,413-422  (1999). 

30.  T.  M.  Kolb,  J.  Lichy,  and  J.  H.  Newhouse,  “Occult  cancer  in 
women  with  dense  breast:  detection  with  screening  US- 
diagnostic  yield  and  tumor  characteristics,”  Radiology  207, 
191-199  (1998). 


6380  APPLIED  OPTICS  /  Vol.  40,  No.  34  /  1  December  2001 


Design  of  near-infrared  imaging  probe  with  the 
assistance  of  ultrasound  localization 


Quing  Zhu,  Nan  Guang  Chen,  Daqing  Piao,  Puyun  Guo,  and  XiaoHui  Ding 


A  total  of  364  optical  source- detector  pairs  were  deployed  uniformly  over  a  9  cm  X  9  cm  probe  area 
initially,  and  then  the  total  pairs  were  reduced  gradually  to  60  in  experimental  and  simulation  studies. 
For  each  source- detector  configuration,  three-dimensional  (3-D)  images  of  a  1-cm-diameter  absorber  of 
different  contrasts  were  reconstructed  from  the/measurements  made  with  a  frequency-domain  system. 
The  results  have  shown  that  more  than  160  source- detector  pairs  are  needed  to  reconstruct  the  absorp¬ 
tion  coefficient  to  within  60%  of  the  true  yplue  and  appropriate  spatial  and  contrast  resolution.  How¬ 
ever,  the  error  in  target  depth  estimated  from  3-D  images  was  more  than  1  cm  in  all  source- detector 
configurations.  With  the  a  priori  target  depth  information  provided  by  ultrasound,  the  accuracy  of  the 
reconstructed  absorption  coefficient  was  improved  by  15%  and  30%  on  average,  and  the  beam  width  was 
improved  by  24%  and  41%  on  average  for  high-  and  low-contrast  cases,  respectively.  The  speed  of 
reconstruction  was  improved  by  ten  times  on  average.  ©  2001  Optical  Society  of  America 
OC7S  codes;  170.0170,  170.3010,  170.5270,  170.7170,  170.3830. 


1.  Introduction 

Recently,  optical  imaging  techniques  based  on  diffu¬ 
sive  near-infrared  (NIR)  light  have  been  employed  to 
obtain  interior  optical  properties  of  human  tissues.1-8 
Functional  imaging  with  NIR  light  has  the  potential 
to  detect  and  diagnose  diseases  or  cancers  through 
the  determination  of  hemoglobin  concentration,  blood 
02  saturation,  tissue  light  scattering,  water  concen¬ 
tration,  and  the  concentration  and  lifetime  of  exoge¬ 
nous  contrast  agents.  Optical  imaging  requires  that 
an  array  of  sources  and  detectors  be  distributed  di¬ 
rectly  or  coupled  through  optical  fibers  on  a  boundary 
surface.  Measurements  made  at  all  source- detector 
positions  can  be  used  in  tomographic  image  recon¬ 
struction  schemes  to  determine  optical  properties  of 
the  medium.  The  frequently  used  geometric  config¬ 
urations  of  sources  and  detectors  are  ring  arrays4*9-11 
and  planar  arrays. 3*12“15  A  ring  array  consists  of 
multiple  sources  and  detectors  that  can  be  distrib¬ 
uted  uniformly  on  a  ring.  Optical  properties  of  the 
thin  tissue  slice  (two-dimensional  slice)  enclosed  by 
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the  ring  can  be  determined  from  all  measurements. 
A  planar  array  can  be  configured  with  either  trans¬ 
mission  or  reflection  geometries.  In  transmission 
geometry,  multiple  detectors  can  be  deployed  on  a 
planar  array,  and  multiple  sources  or  a  single  source 
can  be  deployed  on  an  opposite  plane  parallel  to  the 
detector  plane.  Optical  properties  of  the  three- 
dimensional  (3-D)  tissue  volume  between  the  source 
and  the  detector  planes  can  be  determined  from  all 
measurements.  In  reflection  geometry,  multiple 
sources  and  detectors  can  be  distributed  on  a  planar 
probe  that  can  be  hand-held.3*15  Optical  properties 
of  the  3-D  tissue  volume  at  slice  depths  below  the 
probe  can  be  determined  from  all  measurements. 
The  reflection  probe  configuration  is  desirable  for  the 
imaging  of  brain  and  breast  tissues. 

Although  many  researchers  in  the  field  have  con¬ 
structed  imaging  probes  using  reflection  geome¬ 
try,2*3*15  to  our  knowledge  the  required  total  number 
of  source- detector  pairs  over  a  given  probe  area 
needed  to  accurately  reconstruct  optical  properties 
and  localized  spatial  and  depth  distributions  has  not 
been  addressed  before.  In  this  paper  we  study  the 
relationship  between  the  total  number  of  source- 
detector  pairs  and  the  reconstructed  imaging  quality 
through  experimental  measurements.  Computer 
simulations  are  performed  to  assist  in  understanding 
the  experimental  results. 

Because  the  target  localization  from  diffusive 
waves  is  difficult,  our  group  and  others  have  intro¬ 
duced  use  of  a  priori  target  location  information  pro- 
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vided  by  ultrasound  to  improve  optical  imaging.15"18 
In  this  paper  we  demonstrate  experimentally  that 
the  accurate  target  depth  information  can  signifi¬ 
cantly  improve  the  accuracy  of  the  reconstructed  ab¬ 
sorption  coefficient  and  the  reconstruction  speed  for 
any  optical  array  configuration. 

The  required  total  number  of  source- detector  pairs 
is  also  related  to  the  image  reconstruction  algorithms 
used.  In  this  paper  we  obtained  experimental  mea¬ 
surements  using  a  frequency-domain  system  with  the 
source  amplitude  modulated  at  140  MHz.  In  simu¬ 
lations,  forward  measurements  were  generated  by 
use  of  the  analytic  solution  of  a  photon  density  wave 
scattered  by  a  spherical  inhomogeneity  embedded  in 
a  semi-infinite  scattering  medium.19  In  both  exper¬ 
iments  and  simulations,  linear  perturbation  theory 
within  the  Born  approximation  was  used  to  relate 
optical  signals  at  the  probe  surface  to  absorption  vari¬ 
ations  in  each  volume  element  within  the  sample. 
The  total  least-squares  (TLS)  method20-22  was  used 
to  formulate  the  inverse  problem.  The  conjugate 
gradient  technique  was  employed  to  iteratively  solve 
the  inverse  problem.  Therefore  the  results  we  ob¬ 
tained  are  directly  relevant  to  the  probe  design  with 
reconstruction  algorithms  based  on  the  linear  pertur¬ 
bation  theory  and  can  be  used  as  a  first-order  approx¬ 
imation  if  high-order  perturbations  are  employed  in 
image  reconstructions. 

This  paper  is  organized  as  follows.  In  Section  2 
we  describe  an  analytic  solution  used  to  generate 
simulated  forward  data,  the  Bom  approximation, 
and  the  TLS  method  for  image  reconstruction.  In 
Section  3  we  discuss  the  probe  geometry,  a  frequency- 
domain  NIR  system  used  to  acquire  the  experimental 
measurements,  and  an  ultrasound  subsystem  used  to 
acquire  the  target  depth  information,  computation 
procedures  used  to  obtain  both  simulated  and  exper¬ 
imental  absorption  images.  In  Sections  4  and  5  we 
report  experimental  results  obtained  from  the  dense 
and  sparse  arrays  with  and  without  a  priori  target 
depth  information.  A  high-contrast  example  is 
given  in  Section  4,  and  a  low-contrast  case  is  given  in 
Section  5.  Simulations  are  performed  to  assist  in 
understanding  the  noise  on  the  image  reconstruction. 
Imaging  parameters  evaluated  are  a  —6-dB  width  of 
the  image  lobe,  the  reconstructed  maximum  value  of 
the  absorption  coefficient  and  its  spatial  location,  and 
the  image  artifact  level.  In  Sections  6  and  7  we 
provide  a  discussion  and  summary,  respectively. 

2.  Basic  Principles 

A.  Forward  Model 

In  our  experiments,  forward  measurements  were 
made  with  a  frequency-domain  system  operating  at  a 
140-MHz  modulation  frequency.  In  computer  sim¬ 
ulations,  forward  measurements  were  generated 
from  an  analytic  solution  of  a  photon  density  wave 
scattered  by  a  spherical  inhomogeneity. 19  When  the 
center  of  the  sphere  coincides  with  the  origin  of 
spherical  coordinates,  the  solution  for  the  scattered 


infinite  medium. 


photon  density  wave  C/sc  outside  the  sphere  at  a  de¬ 
tector  position  r  =  (r,0,<t>)  is  of  the  form 

UJr,v)  =  2  {AjjK*#Btr)  +>;(*0“V)]Y(,m(e,c}>)}> 

l,m 


where  jt  and  nt  are  spherical  Bessel  and  Neuman 
functions,  respectively;  Yi,m(6,4>)  are  the  spherical 
harmonics,  and  kout  =  [(-vp.aout  +  jm)/D°'lt]  is  the 
complex  wave  number  outside  the  sphere,  u  is  the 
angular  modulation  frequency  of  the  light  source, 
lxaout  is  the  absorption  coefficient  outside  the  sphere, 
and  Dout  is  the  photon  diffusion  coefficient  outside  the 
sphere  given  by  Dout  =  l/(3|xs'),  where  p.s'  is  the 
reduced  scattering  coefficient  outside  the  sphere. 
The  coefficients  Alm,  determined  by  the  boundary 
conditions,  are 

A/m  =  -  (jvSk0'11/ i?out)ft  jl)(i0UVs)  Y/,m*  (8S,  <|>5) 
Dmtxji(y)ji'(x)  -  D'*yji'(.y)Mx) 

X  Dmlxh\lY(x)jt(y)  -  D™yh?\x)jay)  ’ 

where  *  =  kouta,  y  =  kma,  rs  =  (r„B„4>,)  is  the  source 
position,  h\1^  are  the  Hankel  functions  of  the  first 
kind,  and  j{  and  h\1)f  are  the  first  derivatives  of  ji  and 
h^.  The  analytic  solution  has  the  important  advan¬ 
tage  in  that  it  is  exact  to  all  orders  of  perturbation 
theory  and  thus  can  represent  accurate  measure¬ 
ments. 

We  generalized  the  above  analytic  solution  to  a 
semi-infinite  geometry  by  using  a  method  of  images 
with  extrapolated  boundary  conditions  (see  Fig.  1)« 
A  type  I  boundary  condition  (zero  light  energy  density 
at  the  extrapolated  boundary)  is  used  to  derive  the 
scattered  wave  U3C’ .  To  calculate  the  U3C'  in  semi¬ 
infinite  geometry,  we  use  rsph+  =  (ro+^o+^o-t-)  anc* 
rSph-  =  (r0-,e0-,<l>o-) t0  represent  the  centers  of  the 
sphere  and  the  image  sphere,  respectively.  The  vec¬ 
tors  r+'  =  r  -  rsph+  =  (r+',0+\4>+')  and  r-'  =  r  ” 
rsph-  “  (r_',e_',<{>-')  are  therefore  pointing  to  the 
detector  position  r  =  (r,0,<}>)  from  the  sphere  and  the 
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image  sphere,  respectively.  The  semi-infinite  solu¬ 
tion  of  Uscf  can  be  approximated  as 

l,m 

+jnAkm*r.')]Ya9*',  <t>+')} 

-  £  {Au-L/rf*^/) 

l,m 

+Mk"*rS)]Yt„(9S,< !»/)} 

+  £  {AiMk°utrJ) 

Ijn 


where  G{r„rd,ta)  is  the  Green  function  and  A|xa(rv)  = 
~  M-0  is  the  medium  absorption  variation.11 
pa  is  the  average  value  of  the  medium  absorption 
coefficient.  By  breaking  the  medium  into  discrete 
voxels,  we  obtain  the  following  linear  equations: 

AT 

£Ac'(rdi,rsi,co)  =  £G(rv],rdi)a>)f7mc(r,,J,r,,„oj)  . 

j 

X  [vAli0(ry.)/5]Arv3.  (9) 

When  =  G(r1!/)rdi,<i))Uinc  (rvj,rai,u>)v£irv3/D,  we  ob¬ 
tain  the  matrix  equation  of  Eq.  (9): 


l,m 

+>/(Aoutr.')]y,m(0.',4,.')}i  (3) 

where 


A‘,m+  =  -UvSIC^D^^k^Y^;,^) 

x[  V^iiyVi'M-rPyjt'i. vMx) 
Dmtxh\iy(x)jt(y)  -  Dinyh\l)(x)jt'(y)  ’ 


(4) 


A/im-  =  -UvSk^/inh^rnY^-,  4,,-) 

X  D^vAyVi'to-iPyji’MMx) 

ir*xh\lY(x)My MWWW'  () 

rs  (rs  A  >4>s  )  ar>d  rs  =  (rs  ,0S  ,4>s—)  are  the 
positions  of  the  source  and  the  image  source,  respec- 
tively. 

The  incident  photon  density  wave  at  the  detector 
position  r  has  the  following  form3: 


Ginc(r*,Oj) 


4ir Z)out 
exp(y&outjr  -  r. 


exp(^°ut|r  -  r/j) 
|r  -  r.+| 


p-  r. 


(6) 


■Hie  total  photon  density  at  detector  r  is  a  super¬ 
position  of  its  incident  (homogeneous)  and  scattered 
(heterogeneous)  waves: 


U{  r,w)  =  Umc(r,uj)  +  Uac'(  r,(o),  (7) 


B.  Born  Approximation  for  Reconstruction 

The  Born  approximation  was  used  to  relate  Ux'  (r,w) 
measured  at  the  probe  surface  to  absorptionVaiia- 
tions  in  each  volume  element  within  the  sample.  In 
the  Born  approximation,  the  scattered  wave  that 
originated  from  a  source  at  rs  and  measured  at  de¬ 
tector  rd  can  be  related  to  the  medium  heterogeneity 
Apa  (r  J  by 


[WWApaW  =  [Gad]  MX  1-  (10) 

The  realistic  constrains  on  Aji0  are  (-a  x  back¬ 
ground  pa)  =  <  A pa  <  1,  where  0  <  a  <1. 

The  above  constrains  ensure  that  the  reconstructed 
absorption  coefficient  |i0  =  background  pa  +  Ap0  is 
positive  and  not  unrealistically  higher  than  unity. 
With  M  measurements  obtained  from  all  possible 
source— detector  pairs  in  the  planar  array,  we  can 
solve  N  unknowns  of  Apa  by  inverting  the  matrix  Eq. 
(10).  In  general,  the  perturbation  in  Eq.  (10)  is  un¬ 
derdetermined  {M  <N)  and  ill-posed. 

When  the  target  depth  is  available  from  ultra¬ 
sound,  we  can  set  Apa  of  a  nontarget  depth  equal  to 
zero.  This  implies  that  all  the  measured  perturba¬ 
tions  were  originated  from  the  particular  depth  that 
contained  the  target.  Because  the  number  of  un¬ 
knowns  was  reduced  significantly,  the  reconstruction 
converged  fast. 


C.  Total  Least-Squares  Solution 

To  solve  the  unknown  optical  properties  of  Eq.  (10), 
several  iterative  algorithms  have  been  used  in  the 
literature  including  the  regularized  least-squares 
method10  and  the  TLS  method.20-21  The  TLS  per¬ 
forms  better  than  the  regularized  least-squares 
method  when  the  measurement  data  are  subject  to 
noise  and  the  linear  operator  V?  contains  errors. 
The  operator  errors  can  result  from  both  the  approx¬ 
imations  used  to  derive  the  linear  model  and  the 
numerical  errors  in  the  computation  of  the  operator. 
We  found  that  the  TLS  method  provides  more  accu¬ 
rate  reconstructed  optical  properties  than  the  regu¬ 
larized  least-squares  method,  so  we  adapted  the  TLS 
method  to  solve  the  inverse  problems.  It  has  been 
shown  by  Golub22  that  the  TLS  minimization  is 
equivalent  to  the  following  minimization  problem: 


ll^~ 
W~+~1 


(11) 


where  X  represents  unknown  optical  properties. 
The  conjugate  gradient  technique  was  employed  to 
iteratively  solve  Eq.  (11). 


UJ(  rd,r„u>) 


j  G(rv,rd,ti))G’inc(r„rs,w) 
X  [vAp.0(rv)/D]drv3, 


3.  Methods 

A.  Probe  Design  and  imaging  Geometry 

There  are  two  basic  requirements  to  guide  the  design 
(8)  of  the  NIR  probe.  First,  all  source—  detector  separa- 
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Fig.  2.  Configuration  of  a  dense  array  with  28  optical  sources  and 
13  detectors  as  well  as  six  ultrasound  transducers.  Large  black 
circles  are  optical  detectors,  gray  circles  are  optical  sources,  and 
small  white  circles  are  ultrasound  transducers.  A  1-cm-diameter 
spherical  target  was  located  at  various  depths  in  simulations  and 
experiments. 


tions  have  to  be  as  large  as  1  cm  so  that  diffusion 
theory  is  a  valid  approximation  for  image  reconstruc¬ 
tion.  Second,  because  the  depth  of  a  photon  path  is 
measured  approximately  one  third  to  one  half  the 
source- detector  separation,  the  distribution  of 
source-detector  distances  should  be  from  approxi¬ 
mately  1  to  10  cm  to  effectively  probe  the  depth  from 
approximately  0.5  to  4  cm.  On  the  basis  of  these 
requirements,  we  deployed  a  total  of  28  sources  and 
13  detectors  over  a  probe  area  of  9  cm  X  9  cm  (see  Fig. 
2).  The  minimum  source- detector  separation  in  the 
configuration  is  1.4  cm  and  the  maximum  is  10.0  cm. 
We  call  this  array  a  filled  or  dense  array  (a  term 
adapted  from  ultrasound  array  design).  The  9  cm  X 
9  cm  X  4  cm  imaging  volume  was  discretized  into 
voxels  of  size  0.5  cm  X  0.5  cm  X  1  cm;  therefore  a  total 
of  four  layers  in  depth  was  obtained.  The  target  was 
a  1-cm-diameter  sphere  located  at  different  locations. 
Because  one  of  the  objectives  of  this  study  was  to 
evaluate  the  target  depth  distribution,  the  centers  of 
the  four  layers  in  depth  were  adapted  to  the  target 
depth.  For  example,  if  the  target  depth  was  z  =  3 
cm,  the  centers  of  the  four  layers  were  chosen  as  1,  2, 
3,  and  4  cm,  respectively. 

The  ultrasound  transducers  shown  in  Fig.  2  were 
deployed  simultaneously  on  the  same  probe.  The 
diameter  of  each  ultrasound  transducer  is  1.5  mm 
and  the  spacing  between  the  transducers  is  4  mm, 
except  the  two  located  closer  to  the  optical  detector  in 
the  middle.  The  spacing  between  these  two  trans¬ 
ducers  is  8  mm.  Because  this  study  requires  accu¬ 
rate  target  location  as  a  reference  to  compare  with 
the  reconstructed  absorption  image  location,  six 
transducers  are  used  to  guide  the  spatial  positioning 
of  a  target.  The  target  is  centered  when  the  two 


Fig.  3.  Schematic  of  a  single- channel  optical  data-acquisition  sys¬ 
tem.  A  140.02-MHz  oscillator  is  used  to  drive  the  laser  diode  (780 
nm)  that  delivers  the  light  to  the  medium  through  the  fiber.  The 
detected  signals  are  amplified  and  mixed  with  signals  from  a  140- 
MHz  oscillator.  The  heterodyned  20-kHz  signals  are  amplified, 
filtered,  and  digitized.  The  signals  from  two  oscillators  are  also 
directly  mixed  to  provide  reference  signals.  The  amplitude  and 
phase  of  the  waveform  received  through  the  medium  are  calculated 
from  signals  measured  through  the  medium  and  the  reference. 
PMT,  photomultiplier  tube. 


middle  ultrasound  transducers  receive  the  strongest 
signals.  The  target  depth  is  determined  from  re¬ 
turned  pulse-echo  signals.  In  this  study  we  do  not 
intend  to  provide  ultrasound  images  of  the  target 
with  such  a  sparse  ultrasound  array,  but  we  demon¬ 
strate  the  feasibility  of  using  a  priori  depth  informa¬ 
tion  to  improve  optical  reconstruction. 

In  regard  to  the  image  voxel  size,  there  is  a  trade¬ 
off  between  the  accurate  estimation  of  the  weight 
matrix  W  and  the  voxel  size.  Because  WtJ  is  a  dis¬ 
crete  approximation  of  the  integral 

[  G(ru,rd,w)  Ginc(rv,r3Iw)(v/Z>)di\,3, 


it  is  more  accurate  when  the  voxel  size  is  smaller. 
However,  the  total  number  of  reconstructed  un¬ 
knowns  will  increase  dramatically  with  the  decreas¬ 
ing  voxel  size.  Because  the  rank  of  the  matrix  W  is 
less  than  or  equal  to  the  total  number  of  measure¬ 
ments  [Eq.  (10)],  a  further  decrease  in  voxel  size  will 
not  add  more  independent  information  to  the  weight 
matrix.  We  found  that  a  0.5  cm  X  0.5  cm  X  1  cm 
voxel  size  is  a  good  compromise.  Therefore  we  used 
this  voxel  size  in  image  reconstructions  reported  in 
this  paper. 

B.  Experimental  System 

We  constructed  a  NIR  frequency-domain  system,  and 
the  block  diagram  of  the  system  is  shown  in  Fig.  3.  On 
the  source  side,  a  140.000-MHz  sine-wave  oscillator 
was  used  to  modulate  the  output  of  a  780-nm  diode 
laser  that  was  housed  in  an  optical  coupler  (OZ  Op¬ 
tics  Inc.).  The  output  of  the  diode  was  coupled  to  the 
turbid  medium  through  a  single  200-|xm  multimode 
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optic  fiber  On  the  reception  side,  an  optical  fiber  of 
4  mm  m  diameter  was  used  to  couple  the  detected 
'  fTu  t0  f  photomultlPher  tube  detector.  The  output 
.  . tha  Photomultiplier  tube  was  amplified  and  then 

Tun-r  a,  °Ca!  oscihator  at  a  frequency  of  140.020 
MHz.  The  heterodyned  signal  at  20  kHz  after  the 
mixer  was  further  amplified  and  filtered  by  a  band- 
P3jS  ,  The  outPuts  of  two  oscillators  (140.000- 
and,  140-020-MHz  signals)  were  directly  mixed  to 
produce  20-kHz  reference  signals.  Both  signals 
were  sampled  simultaneously  by  a  dual-channel  1.25- 
MHz  analog-to-digital  converter  (A/D)  board.  The 
Hilbert  transform  was  performed  on  both  sampled 
and  reference  waveforms.  The  amplitude  of  the  Hil¬ 
bert  transform  of  the  sampled  waveform  corresponds 
to  the  measured  amplitude,  and  the  phase  difference 
between  the  phases  of  the  Hilbert  transforms  of  the 
sampled  and  reference  waveforms  corresponds  to  the 
measured  phase. 

A  black  probe  with  holes  shown  in  Fig.  2  was  used 
to  emulate  the  semi-infinite  boundary  condition. 

two /-I)  positioners  were  moved  independently  to 

position  the  source  and  detector  fibers  at  the  desired 
spatial  locations  within  the  9  cm  x  9  cm  area 
A  challenge  in  the  reflection  NIR  probe  design  is  to 
preserve  a  huge  dynamic  range  in  received  signals 
For  example,  the  amplitude  at  a  1-cm  source- 
detector  separation  measured  from  0.6%  Intralipid  in 
reflection  mode  is  approximately  84  dB  larger  than 
that  at  a  9-cm  separation.  So  the  signals  can  be 
saturated  when  they  are  measured  from  closer 
source- detector  pairs,  but  they  may  be  too  low  at 
more  distant  source-detector  pairs.  The  problem 
can  be  overcome  by  means  of  controlling  the  light 
llluminatmn.  At  least  two  illumination  conditions 
need  to  be  used:  a  low  source  level  for  closer  source- 
detector  pairs  and  a  high  level  for  distance  source- 
detector  pairs  In  our  system,  a  30-dB  attenuator 
connected  to  the  140-MHz  oscillator  was  switched  on 
and  off  to  provide  two  different  source  levels  and  thus 
to  preserve  the  dynamic  range.  Figure  4(a)  shows  a 
plot  of  the  measured  log  [p‘£/sd(p)]  versus  the  source- 
detector  separation  p,  and  Fig.  4(b)  shows  the  plot  of 
the  measured  phase  versus  the  source- detector 
separation.  The  Intralipid  concentration  was 
U.6%,  which  corresponded  to  p./  =  6  cm-1.  Theo- 
retically  both  log  [p2Dsd(p)]  and  phase  are  linearly 
related  to  the  source-detector  separation  because 
of  the  semi-infinite  boundary  condition  used,3  and 
experimental  measurements  shown  in  Fig.  4  vali¬ 
date  that  they  are  linearly  related  to  the  source- 
detector  separation. 

The  ultrasound  system  consists  of  six  transducers 
(see  big.  5),  a  pulser  (Panametrics  Inc.),  an  A/D  con¬ 
verter,  and  a  multiplexer.  The  pulser  provided  a 
high-voltage  pulse  of  6-MHz  central  frequency  to  drive 
each  selected  transducer.  The  returned  signals  were 
received  by  the  same  transducer,  amplified  by  the  re- 
ceiving  circuit  inside  the  pulser,  and  sampled  by  the 
A/D  converter  with  a  100-MHz  sampling  frequency. 


Fig.  4.  Calibration  cunres.  (a)  log  [p’U(p)]  versus  source- 
detector  separation,  (b)  Phase  versus  source-detector  separa- 


C.  Computation  Procedures 

1.  Computation  Procedures  of  Experimental  Data 
To  study  the  relationship  between  the  total  number 
of  source- detector  pairs  and  distributions  of  recon¬ 
structed  optica]  absorption  coefficients,  we  started 
from  the  dense  array  with  a  total  of  28  sources  and  13 
detectors  (see  Fig.  2)  and  gradually  reduced  this 
number  to  generate  sparse  arrays  with  24  X  13  (24 
sources  and  13  detectors),  20  X  13,  28  X  9  24  X  9 

91620  X  62°ffiXx%12  Xi?o16  X  9’  28  X  5’  24  *  5- 12  * 

y,  AO  x  5  16  X  5,  and  12  X  5  source- detector  pairs 
respectively.  Each  sparse  array  was  a  subset  of  the 
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multiplexer 
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Fig.  5.  Ultrasound  data-acquisition  system.  The  pulser  is  used 
to  generate  high-voltage  pulses  that  are  used  to  excite  the  selected 
ultrasound  transducer.  The  returned  signals  are  received  by  the 
selected  transducer  and  are  sampled  by  the  A/D  converter. 
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dense  array,  and  its  probe  area  was  the  same  as  the 
dense  array!  For  each  sparse  array  configura  1  , 

we  compared  the  reconstructed  optical  ™a.f?n^Pa 
rameters  measured  from  the  dense  array  wit  os 
from  the  sparse  arrays.  The  parameters  include  the 
maximum  values  of  reconstructed  jia  at  different  lay¬ 
ers  and  their  spatial  locations,  spatial  resolution  and 
artifact  level  of  the  |ia  image,  and  target  depth  dis¬ 
tribution.  Targets  of  different  absorption  contrasts 
were  located  at  different  positions  inside  the  In¬ 
tralipid.  For  each  target  case,  one  set  of  measure¬ 
ments  with  the  dense  array  was  obtained,  and 
subsets  of  the  measurements  were  used  as  sparse 
array  measurements.  In  all  experiments,  the  back¬ 
ground  Intralipid  concentration  was  approximately 
0.6%,  and  iV  was  experimentally  determined  from 
curve  fitting  results.  Currently,  we  did  not  recon¬ 
struct  target  fi/,  and  we  used  the  common  for 
both  the  background  and  the  target. 

The  total  number  of  iterations  or  stopping  criterion 
was  difficult  to  determine  for  experimental  data. 
Ideally,  the  iteration  should  stop  when  the  object 
function  [see  Eq.  (11)]  or  the  error  performance  sur¬ 
face  reaches  the  noise  floor.  However,  the  system 
noise,  particularly  coherent  noise,  was  difficult  to  es¬ 
timate.  In  general,  we  found  that  the  reconstructed 
values  were  closer  to  true  values  when  the  object 
function  reached  approximately  5-15%  of  the  initial 
value  (total  energy  in  the  measurements).  How¬ 
ever,  this  criterion  was  applicable  only  to  reconstruc¬ 
tions  with  total  source- detector  pairs  closer  to  the 
dense  array  case.  Therefore  we  used  this  criterion 
for  the  dense  array  reconstruction  and  used  the  same 
iteration  number  obtained  from  the  dense  array  for 
the  sparse  arrays.  Thus  the  iteration  number  is 
normalized  to  the  dense  array  case. 


2.  Computation  Procedures  of  Simulation 

Simulations  were  performed  to  assist  the  under¬ 
standing  of  the  random  noise  on  the  reconstructed 
absorption  coefficient.  In  simulations,  Gaussian 
noise  with  different  standard  deviations  proportional 
to  the  average  value  of  each  forward  data  set  was 
added  to  the  forward  measurements.  Typically, 
0.5%,  1.5%,  and  2.0%  of  the  average  value  of  each 
forward  data  set  were  used  as  standard  deviations  to 
generate  noise.  In  simulations,  the  target  p.a  was 
changed  to  different  contrast  values,  and  target  \x 
was  kept  the  same  as  the  background.  The  back¬ 
ground  |xn  and  \is'  were  0.02  and  6  cm  ,  respectively. 

The  stopping  criterion  used  in  the  simulations  was 
based  on  the  noise  level  of  the  object  function.  Con¬ 
sidering  that  the  object  function  fluctuates  within  one 
standard  deviation  a  around  the  mean  E  when  the 
iteration  number  is  large,  we  can  use  E  +  or  as  a 
stopping  criterion,  i.e.,  the  iteration  will  stop  if  the 
object  function  is  less  than  E  +  cr.  When  the  linear 
perturbation  is  assumed,  E  can  be  approximated  as 
Xfn(j)2  and  a  as  2?[2n(j)4]1/2,  where  N  is  the  total 
number  of  source- detector  pairs  and  n(j)  is  the  gen¬ 
erated  random  noise  with  a  standard  deviation  pro¬ 


portional  to  the  specified  percentage  of  the  mean  of 
the  forward  data  set. 


D.  Testing  Targets 

Spherical  testing  targets  of  ~1  cm  in  diameter  were 
made  of  acrylamide  gel.16  The  acrylamide  powder 
was  dissolved  in  distilled  water,  and  20%  concentra¬ 
tion  of  Intralipid  was  added  to  the  acrylamide  solu¬ 
tion  to  dilute  the  solution  to  a  0.6%  Intralipid 
concentration  (m'  =  6  cm"1).  India  ink  was  added 
to  the  solution  to  produce  target  pa  of  different  val¬ 
ues.  Acoustic  scattering  particles  of  200  p-m  in  di¬ 
ameter  were  added  to  the  solution  before 
polymerization.  Components  of  ammonium  persul¬ 
fate  and  tetramethylethylenediamine  (known  as  TE- 
MED)  were  added  to  the  solution  to  produce 
polymerization. 


4.  Results  of  a  High-Contrast  Target  Case 

A.  Experimental  Results  of  a  Dense  Array 
Figure  6(a)  is  an  experimental  image  of  a  high- 
contrast  target  (m-q  =  0.25  cm"1)  located  in  the  In¬ 
tralipid  background  (p/  =  6  cm  x).  The  target  was 
a  1-cm-diameter  sphere  and  its  center  was  located  at 
(*  =  o,  y  =  0,  z  =  3.0  cm),  where  x  and  y  were  the 
spatial’ coordinates  and  z  was  the  propagation  depth. 
The  target  depth  was  well  controlled  by  use  of  ultra¬ 
sound  pulse-echo  signals,  and  the  error  was  less  than 
1  mm.  The  3-D  images  were  reconstructed  from  the 
measurements  made  with  the  dense  array,  and  the 
image  shown  was  obtained  at  target  layer  3.  The 
centers  of  the  imaging  voxels  in  z  are  1,  2,  3,  and  4  cm 
for  layers  1,  2,  3,  and  4,  respectively.  The  measured 
maximum  value  of  the  image  lobe  [|ia(max)3 was  °-233 
which  was  a  close  estimate  of  the  target  pLa. 
The  measured  spatial  location  of  |ia(max) was  (*  =  ?  3 
crxij  y  =  o.O  cm),  which  agreed  reasonably  well  with 
the  true  target  location.  The  spatial  resolution  can 
be  estimated  from  the  -6-dB  contour  plot  of  Fig.  6(a) 
which  is  shown  in  Fig.  6(b).  The  outer  contour  is  —6 
dB  from  the  jlG(max),  and  the  contour  spacing  is  1  dB^ 
The  width  of  the  image  lobe  measured  at  the  -b 
dB-level  corresponds  to  a  full  width  at  half-maximum 
(FWHM),  which  is  commonly- used  to  estimate  reso¬ 
lution.  The  measured  widths  of  longer  and  shorter 
axes  were  1.01  and  1.60  cm,  respectively,  and  t  e 
geometric  mean  was  1.27  cm,  which  was  used  to  rep¬ 
resent  the  -6-dB  beam  width.  The  contrast  resolu¬ 
tion  can  be  estimated  from  the  peak  artifact  level, 
and  no  artifact  was  observed  in  the  image.  The  tar¬ 
get  depth  can  be  assessed  from  the  images  obtained 
from  other  nontarget  layers.  Figure  6(c)  is  the  im¬ 
age  obtained  at  nontarget  layer  4,  and  an  image  lobe 
of  ll  ,  ,  =  0.138  cm"1  was  observed.  The  spatial 

locatioiT of  Patmaxi  was  (x  =  0.0,  y  =  0.0),  wluch 
agreed  well  with  the  true  target  location.  No  dis¬ 
tinct  lobes  were  observed  at  nontarget  layers  1  and  2, 
which  indicates  that  the  error  in  the  target  depth 
estimated  from  3-D  images  was  approximately  1cm. 
Because  the  error  in  the  true  target  depth  was  less 
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(b) 


nf  9«6  v  f,XP“al  3'D!mageS  of  reconstructed  with  a  total 
ol  28  x  13  =  364  source-detector  pairs  at  2437  iterations.  The 
terget  (jx„  -  0.25  cm  )  was  located  at  (x  =  0,  y  =  0,  z  =  3  0  cm) 
inside  the  Litrahpid  background,  (a)  Reconstructed  4  at  target 
layer  3.  The  horizontal  axes  represent  spatial  x  andy  coordinates 
m  centimeters,  and  the  vertical  axis  is  the  |i„.  The  measured 
maximum  value  of  the  image  lobe  [|ia(m„,]  was  0.233  cm~\  and  its 
ocation  was  (x  =  0.5, y  =  0.0).  No  image  artifacts  were  observed, 
(b)  6-dB  contour  plot  of  (a).  The  outer  contour  is  -6  dB  from  the 
and  the  contour  spacing  is  1  dB.  The  measured  -6-dB 
beam  width  was  1.27  cm.  (c)  Reconstructed  ft,  at  nontarget  layer 

t  °f3tr,ength  0138  cm"‘  and  spatial  location  of 

v x  —  0.0,  y  —  0.0)  was  observed. 

than  1  mm,  this  1-cm  error  was  due  largely  to  the 
depth  uncertainty  of  diffusive  waves. 

B.  Simulation  Results  of  a  Dense  Array 

Our  simulations  support  the  experimental  results 
figure  7  shows  simulation  results  obtained  with  the 
ense  array.  A  simulated  1-cm-diameter  absorber 
'•P'«  7  .  5  cni  )  was  located  at  (x  =  0,  y  =  0  z  =  3 
cm)  inside  a  homogeneous  scattering  background  (u. ' 

,  6  c“  )'  added  0.5%  Gaussian  noise  to  the 

TnZ?/d  Jfta  gjnTated  from  the  analytic  solution. 
Images  obtained  at  nontarget  layer  2,  target  layer  3 
and  nontarget  layer  4  are  shown  in  Fig.  7(a)  Fie.  (b)’ 
and  7(c),  respectively.  No  target  was  found  at  layer 


(d> 

noSy\'n  5‘““lated  3*®  lma6es  of  |i„  reconstructed  with  a  total  of 
28  x  13  -  364  source-detector  pairs.  The  target^  =  0.25  cm-1! 
was  located  at  (x  =  0,  y  =  0,  z  =  3.0  cm,  inside  the  Intralipid 
background,  (a)  Reconstructed  |i0  at  nontarget  layer  2  (simula¬ 
tion,  0.5%  noise).  No  image  lobe  was  observed,  (b)  Recon- 
s  ructed  |i„  at  target  layer  3  (simulation,  0.5%  noise).  The  target 
o  s  rengt  M.a(mar)  =  0.248  cm  1  and  the  spatial  location  (0  0  0  0> 

tT  °nf*red'  !C>  ?feonstructed  ^  at  nontarget  layer  4  (simula¬ 
tion,  0.5%  noise).  The  target  of  strength  0. 190  cm  * 1  and  location 
x  0.0,  y  -  0.0)  was  observed,  (d).  Reconstructed  d,  at  non- 
target  layer  2  with  1.0%  noise.  The  target  o [  0.028  cm“l  was 
observed.  Note  that  the  scale  of  (d)  is  different  from  fai-(c). 

2,  and  the  target  of  strengths  0.248  and  0.190  cnT1 
appeared  at  layers  3  and  4,  respectively.  However 
when  the  noise  level  in  the  forward  data  was  in¬ 
creased  to  1.0%,  the  target  of  strength  {Lalmax)  =  0.028 
cm  appeared  at  nontarget  layer  2  [Fig.  7(d)]  as  well 
as  target  layer  3  [|ia(max)  =  0.163  cm"1] and  nontarget 
layer  4  [p.o(max)  -  0.111  cm  x].  This  suggests  that 
measurement  noise  is  an  important  parameter  to  af- 
tect  the  target  depth  estimate. 

C.  Experimental  Results  of  a  Dense  Array  with 
Ultrasound  Assistance 

From  ultrasound  we  obtained  the  target  depth  as 
well  as  the  target  boundary  information.  Figure 

H?n?hr o  o  receiTed  Pulse-echo  signals  (from  a 

depth  of  2.3-3.5  cm)  obtained  from  the  six  ultrasound 
ransducers  (see  Fig.  2).  The  spatial  dimension  cov¬ 
ered  by  the  six  transducers  is  2.4  cm.  The  front 
surface  of  the  target  was  indicated  clearly  by  the 
returned  pulses  shown  by  two  arrows,  and  the  back 
surface  wasseen  through  the  reflection  of  a  soft  plas¬ 
tic  plate.  The  reflected  signals  of  the  back  surface 
are  also  shown  by  arrows.  The  plate  was  used  to 
hold  the  target  and  was  transparent  to  light.  The 
measured  depth  of  the  target  front  surface  was  2  44 
cm  and  the  back  surface  was  3.43  cm.  The  distance 
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Fig.  8.  (a)  Ultrasound  pulse-echo  signals  or  A-scan  lines  obtained 
from  six  transducers.  The  abscissa  is  the  propagation  depth  in 
millimeters.  From  reflected  signals,  the  measured  depth  of  the 
target  front  surface  is  2.44  cm,  and  the  back  surface  is  3.43  cm. 
The  center  of  the  target  is  -3  cm.  The  total  length  of  the  signal 
corresponds  to  1.2  cm  in  depth,  and  the  measured  distance  be¬ 
tween  the  front  and  the  back  surfaces  is  0.993  cm.  The  spatial 
dimension  covered  by  the  transducers  is  2.4  cm.  (b)  An  image  of 
the  high-contrast  target  reconstructed  at  a  target  layer  when  we 
used  only  a  priori  target  depth  information  provided  by  ultra¬ 
sound.  The  reconstructed  |io{max)  reached  0.245  cm-1  at  216  it¬ 
erations. 

between  the  peaks  of  front  reflection  and  backreflec- 
tion  was  0,993  crn,  which  corresponded  to  the  target 
size.  With  the  assistance  of  target  depth,  we  recon¬ 
structed  the  absorption  image  at  the  target  layer 
only.  Figure  8(b)  is  the  reconstructed  jl0  obtained 
from  the  dense  array  measurements.  A  total  of  216 
iterations  were  used  to  obtain  (l^max)  ~  0.245  cm-1, 
and  the  reconstruction  was  approximately  ten  times 
faster  than  that  without  the  depth  information.  The 
spatial  resolution  was  approximately  the  same  as 
Fig.  6(a),  and  the  -6-dB  beam  width  was  1.31  cm. 
The  contrast  resolution  was  5  dB  poorer  because  the 
measurement  noise  was  lumped  to  single-layer  re¬ 
construction  instead  of  being  distributed  to  four  lay¬ 
ers. 

D.  Experimental  Results  of  a  Sparse  Array 
The  imaging  quality  of  sparse  arrays  decreased.  Fig¬ 
ure  9(a)  is  an  experimental  image  at  target  layer  3 
obtained  from  the  16  X  5  sparse  array.  The  sparse 
array  measurements  used  were  a  subset  of  the  dense 


Fig.  9.  Experimental  images  at  target  layer  3  reconstructed  with 
a  total  of  16  X  5  =  80  source-detector  pairs,  (a)  Reconstructed  fia 
at  target  layer  3  with  2478  iterations.  The  measured  jla(max)  was 
0.107  cm-1,  which  was  43%  of  the  true  value,  and  the  spatial 
location  of  £a(max)  was  (r  =  0.5,  y  =  -0.5).  The  measured  peak 
image  artifact  level  was  -10  dB  below  the  peak  of  the  main  image 
lobe,  (b)  -12-dB  contour  plot  of  (a).  The  outer  contour  is  -12 
dB,  and  the  contour  spacing  is  2  dB.  The  measured  -6-dB  beam 
width  was  2.55  cm,  which  was  200%  broader  than  that  of  the  dense 
array,  (c)  Reconstructed  |ia  at  target  layer  3  with  10,000  itera¬ 
tions.  The  measured  |ia(max)  reached  0.264  cm"1,  and  the  peak 
artifact  level  was  increased  by  2  dB  as  well,  (d)  Reconstructed  jla 
at  the  target  layer  with  only  a  priori  target  depth  information 
provided  by  ultrasound.  jxa(max)  —  0.173  cm"1  at  216  iterations. 


array  measurements.  The  measured  jia(max)  was 
0.107  cm"1,  which  was  43%  of  the  true  value.  The 
—6-dB  beam  width  was  2.55  cm,  which  was  200% 
broader  than  that  of  the  dense  array.  Edge  artifacts 
were  observed  and  are  best  seen  from  Fig.  9(b),  which 
is  the  -12-dB  contour  plot  of  Fig.  9(a).  The  outer 
contour  is  - 12  dB,  and  the  spacing  is  2  dB.  The  peak 
artifact  level  is  -10  dB  from  the  |ia(max).  The  recon¬ 
structed  |ia  can  be  increased  if  the  iteration  is  signifi¬ 
cantly  increased.  The  iteration  number  used  to 
obtain  Fig.  9(a)  was  2478,  which  was  the  same  as  that 
used  to  obtain  Fig.  6.  When  the  iteration  number  was 
increased  to  10,000  for  the  sparse  array  case,  the  re¬ 
constructed  |i a(max)  reached  0.264  cm"1,  which  was 
close  to  the  true  target  fia.  However,  the  image  arti¬ 
fact  level  was  increased  too  [see  Fig.  9(c)],  and  the  ratio 
of  the  peak  artifact  to  |ia(max)  was  2  dB  higher  than 
that  shown  in  Fig.  9(a).  In  addition,  the  background 
noise  fluctuation  of  nontarget  layers  1  and  2  was  in¬ 
creased.  The  noise  fluctuation  can  be  estimated  from 
the  standard  deviations  of  reconstructed  |ia  at  nontar- 
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get  layers  1  and  2.  At  2478  iterations,  the  averages 
and  the  standard  deviations  of  (i0  were  0.0233 
(±0.0028)  and  0.0202  cm-1  (±0.0022  cm-1)  for  non¬ 
target  layers  1  and  2,  respectively,  and  these  values 
were  0.0244  (±0.0078)  and  0.0208  cm-1  (±0.0073 
cm-1)  at  10,000  iterations. 

The  measured  maximum  values  of  the  image  lobes 
or  target  strengths  at  the  target  layer  and  nontarget 
layer  4  continuously  grew  with  each  iteration,  even 
though  the  reconstructed  values  at  the  target  layer 
were  close  to  the  true  value.  This  problem  was  men¬ 
tioned  in  the  literature,11  but  was  not  explained  well. 
It  is  largely  related  to  use  of  nonlinear  constrains  on 
Ap,a  in  Eq.  (10),  particularly  the  choice  of  a.  When  a 
is  close  to  1,  the  reconstruction  converges  fast,  and 
the  target  strength  increases  little  after  a  certain 
number  of  iterations.  When  a  is  close  to  0,  the  re¬ 
construction  converges  slowly,  and  the  target 
strength  grows  continuously.  However,  when  the 
measurement  signal-to-noise  ratio  (SNR)  is  not  high, 
for  example,  in  sparse  array  or  low-contrast  cases, 
the  choice  of  a  ***  1  can  cause  unstable  reconstruction. 
In  some  cases,  the  reconstructed  images  can  jump 
from  one  set  of  fiQ  to  another,  which  causes  the  object 
function  to  increase  suddenly  and  reduce  again. 
This  is  related  to  the  underdetermined  nature  of  Eq. 
(10),  i.e.,  the  unknowns  are  far  more  than  the  mea¬ 
surements.  In  some  cases  the  reconstructed  images 
have  multiple  lobes  of  similar  strengths,  which  indi¬ 
cate  that  the  reconstruction  does  riot  converge  at  all. 
In  all  cases,  the  image  background  fluctuations  were 
large  compared  with  the  fluctuations  when  a  smaller 
a.  was  used.  We  found  that  a  between  0.1  and  0.4 
can  provide  stable  reconstruction,  and  we  used  a  «= 
0.1  for  all  experiments. 

Another  factor  that  accounts  for  the  slow  increase 
in  the  reconstructed  value  is  use  of  linear  perturba¬ 
tion  to  approximate  the  measurements  that  contain 
all  higher-order  perturbations.  The  minimization 
procedure  [Eq.  (11)]  blindly  minimizes  the  difference 
between  the  measurements  and  their  linear  approx¬ 
imation  WAp,a  and  therefore  reconstructs  higher  and 
higher  A|x0  if  the  iteration  continues. 

The  target  depth  estimate  was  poorer  than  that  of 
the  dense  array  because  of  the  lower  SNR  of  the 
sparse  array  measurements.  Similar  to  the  dense 
array  case,  a  target  of  strength  |ia(max)  =  0.072  cm"1 
and  location  ( x  —  0.5,  y  =  0.0)  was  observed  at  non¬ 
target  layer  4.  In  addition,  a  target  of  strength 
M'o(max)  —  0.0348  cm  *  and  location  (x  —  0.5,  y  = 
-0.5)  was  observed  at  nontarget  layer  2.  However, 
the  target  mass,  which  was  approximately  the  vol¬ 
ume  underneath  the  image  lobe,  was  much  smaller 
than  that  obtained  at  layers  3  and  4,  and  the  target  at 
layer  2  was  buried  in  the  background  noise. 

Figure  9(d)  is  the  reconstructed  |ia  at  the  target 
layer  from  only  the  sparse  array  measurements.  A 
total  of  216  iteration  steps  were  used  to  obtain 
t*,a(max)  —  0.173  cm  *,  and  the  reconstruction  was 
approximately  50  times  faster  than  that  without  the 
depth  information.  The  measured  -6-dB  beam 
width  was  1.49  cm,  which  was  approximately  the 
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Fig.  10.  |ia(max)  versus  the  total  number  of  source— detector  pairs. 
The  center  of  the  target  (|ia  =  0.25  cm”1)  was  located  at  (jc  =  0.0, 
y  ~  0.0,  *  «  3.0  cm)  in  computer  simulations  and  experiments,  (a) 
Curves  were  obtained  at  the  target  layer.  Two  dashed  curves 
(upper  and  lower)  are  the  curve-fitting  results  of  simulation  data 
points  obtained  with  0.5%  and  2.0%  noise  added  to  the  forward 
data,  respectively.  The  experimental  data  are  plotted  with  cir¬ 
cles,  and  the  dashed  curve  in  the  middle  is  the  fitting  result  of  the 
experimental  points,  (b)  The  measured  target  strength  (circles) 
and  the  curve-fitting  result  (lower  curve).  The  measured  target 
strength  (stars)  was  reconstructed  at  the  target  layer  only  by  use 
of  a  priori  depth  information  and  the  curve  fitting  result  (upper 
curve). 


same  as  Fig.  9(c)  but  was  improved  42%  from  Fig, 
9(a).  The  contrast  resolution  was  2  dB  worse  for  the 
same  reason  discussed  above. 

E.  Simulation  and  Experimental  Results  of  Reconstructed 
f^a(max)  versus  Total  Number  of  Source— Detector  Pairs 

To  understand  the  effects  of  random  noise  on  the 
performance  of  the  reconstruction,  we  performed 
simulations  for  each  array  configuration.  Gauss¬ 
ian  noise  of  0.5%,  1.0%,  1.5%,  and  2.0%  were  added 
to  each  fonvard  data  set  generated  from  the  ana¬ 
lytic  solution  [see  Eq.  (3)].  The  center  of  a  simu¬ 
lated  1-cm-diameter  spherical  target  (|ia  =  0.25 


cm  1)  was  located  at  (x  -  0,  y  —  0,  z  -  3.0  cm). 
Reconstructed  images  at  different  noise  levels  were 
obtained  for  each  array  configuration,  and  the  peak 
values  of  image  lobes  [|ia(maX)]  at  target  layer  3  were 
measured.  Figure  10(a)  shows  simulation  and  ex¬ 
perimental  results  of  reconstructed  |ia(m^)  versus 
the  total  number  of  source- detector  pairs.  Two 
dashed  curves  are  the  fitting  results  of  simulated 
data  points  with  0.5%  (upper)  and  2.0%  (lower) 
noise.  The  dashed  curve  in  the  middle  is  the  fitting 
result  of  experimental  points  plotted  with  circles. 
Second-order  polynomials  were  used  for  all  curve 
fittings.  The  reduction  of  the  reconstructed  |ia(max) 
was  significant  when  the  noise  level  went  up  in  the 
simulated  data.  Because  the  SNR  of  experimental 
data  was  decreased  when  the  total  number  of 
source- detector  pairs  was  reduced,  the  data  were 
scattered  around  the  0.5%  noise  curve  when  the 
total  pairs  were  large  and  were  distributed  around 
the  2.0%  noise  curve  (values  were  60%  less  than 
that  obtained  from  the  dense  array)  when  the  total 
pairs  were  reduced  to  less  than  140.  However,  be¬ 
cause  our  experimental  system  has  both  coherent 
and  random  noise,  simulations  based  on  random 
noise  can  only  qualitatively  explain  the  noise  effect 
on  the  experimental  data. 

Figure  10(b)  shows  the  experimental  results  of  re¬ 
constructed  jxa(max)  versus  the  total  number  of 
source— detector  pairs  obtained  from  3-D  imaging 
(circles)  and  ultrasound-assisted  imaging  (stars). 
The  upper  curve  is  the  fitting  result  of  the  circles,  and 
the  lower  curve  is  the  result  of  the  stars.  In  both 
cases,  the  reconstructed  values  were  decreased  when 
the  total  number  of  pairs  was  reduced.  However, 
the  reconstructed  values  were  more  accurate  when 
the  target  depth  information  was  available,  and  the 
improvement  on  average  was  15%.  The  improve¬ 
ment  was  more  dramatic  when  the  total  pairs  were 
less. 

For  sparse  arrays  with  total  source- detector 
pairs  less  than  140,  the  reconstructed  |ia(max)  could 
be  increased  if  the  iterations  were  significantly  in¬ 
creased.  However,,  the  image  artifact  level  of  the 
target  layer  and  the  noise  level  of  the  nontarget 
layers  were  increased  too.  With  the  assistance  of 
simulations,  we  offer  the  following  explanations  to 
the  increased  image  artifact  and  noise  level  prob¬ 
lem.  In  simulations,  the  iteration  was  stopped 
when  the  TLS  error  between  the  measurement  and 
the  linear  approximation  reached  the  noise  floor, 
which  was 

E  +  c  =  i  n{jf  +  £  [2 nurr2, 

j  J 

where  n(j)  was  the  generated  random  noise  with  the 
standard  deviation  proportional  to  a  certain  percent 
of  the  mean  of  the  forward  data  set  for  each  array 


configuration.  When  the  object  function  reached  the 
noise  floor,  the  gradient 

-2(?73d-WA|xa)r(W) 

8  Ap.aT  A|xa  +  1 

-2(£7,  -  ~  WAivHApJ 

(Ap./  A(ia  +  l)2 

can  be  approximated  as 

-2(iV)r(W)  —  2  (jV)T(iST)  ( A  (j.a) 

V' 8  A|xarA|xa  +  1  (Ap./A|j.a  +  l)2  ’ 

where  N  is  the  noise  vector.  Therefore  the  search 
procedure  of  Eq.  (11)  is  more  random  and  noisy.  Be¬ 
cause  =  Ana<okv  +  pVg,  the  Ajia  updating  is 

more  random  and  noisy.  (3  is  proportional  to  the 
square  of  the  gradient.  Continuous  iteration  when 
the  object  function  has  reached  the  noise  floor  may 
destroy  the  convergence  of  the  reconstruction.  We 
found  that  when  the  SNR  of  the  data  is  high,  for 
example,  high-contrast  cases,  continuous  iteration  in 
general  increases  reconstructed  p,a  and  sidelobes. 
However,  when  the  SNR  of  the  data  is  low,  for  exam¬ 
ple,  low-contrast  cases,  continuous  iteration  does  not 
increase  the  reconstructed  but  destroys  the  con¬ 
vergence  of  the  reconstruction  [see  Fig.  12(c)  below]. 

F.  Experimental  Results  of  Imaging  Parameters  Versus 
Total  Number  of  Source-Detector  Pairs 
The  imaging  parameters  measured  from  different  ar¬ 
ray  configurations  are  listed  in  Table  1.  Listed  first 
are  the  measured  parameters  at  target  layer  3. 
These  parameters  are  a  —  6-dB  beam  width  of  the 
image  lobe,  the  peak  sidelobe  level,  |io(max)>  and  the 
distance  in  the  x— y  plane  between  the  location  of 
£a(max)  and  the  true  target  location.  Next  are  the 
same  parameters  measured  at  nontarget  layer  4. 
The  increase  in  beam  width  was  negligible  for  the 
arrays  with  more  than  140  source— detector  pairs  and 
was  100%  broader  for  the  sparse  arrays  with  total 
pairs  less  than  this  number.  The  sidelobe  level  was 
progressively  increased  when  the  total  pairs  were 
reduced.  At  nontarget  layer  4,  when  the  total  pairs 
were  reduced  to  less  than  140,  the  measured  image 
lobes  were  broad  and  no  sidelobes  were  seen.  The 
agreement  between  the  measured  |la(max)  location 
and  the  true  target  location  is  good  for  all  the  array 
configurations,  which  suggests  that  this  parameter  is 
not  sensitive  to  the  total  number  of  source- detector 
pairs  in  high-contrast  cases.  Table  1  next  lists  the 
strength  of  the  target  measured  at  nontarget  layer  2 
and  its  spatial  location.  Because  the  SNR  of  the 
sparse  array  measurement  was  lower,  the  target  ap¬ 
peared  at  layer  2  when  total  pairs  were  less  than  180. 
However,  in  all  cases,  the  target  mass  measured  at 
this  layer  was  much  smaller  than  that  obtained  at  the 
target  layer  and  nontarget  layer  4.  Finally,  Table  1 
lists  the  measured  imaging  parameters  when  the  tar¬ 
get  depth  was  available  to  optical  reconstruction. 
Compared  with  parameters  obtained  from  optical  im- 
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Fig.  11.  Experimental  images  of  ji0  reconstructed  from  a  total  of 
28  X  13  =  364  source- detector  pairs.  The  target  (p.a  =  0.10  cm  ) 
was  located  at  (x  =  0,  y  =  0,  z  =  2.5  cm)  inside  the  Intralipid 
background,  (a)  Reconstructed  |ia  at  target  layer  3.  The  mea¬ 
sured  |ia(max)  was  0.063  cm-1  at  510  iterations,  and  its  location  was 
(x  —  0.0,  y  -  —0.5).  Edge  artifacts  were  observed  and  the  peak 
level  was  -7  dB  from  the  <ia{max).  (b)  Reconstructed  at  non- 
target  layer  4.  The  measured  |lfl{max)  was  0.0871  cm  at  510 
iterations,  and  its  location  was  (x  -  0,  y  =  0).  (c)  Reconstructe 

(ia  at  the  target  layer  when  only  the  target  depth  information 
provided  by  ultrasound  was  used.  p.a(max)  =  0107  cm  at  56 
iterations. 

aging  only,  the  -6-dB  beam  width  was  improved  by 
24%  on  average,  and  the  reconstruction  speed  was 
approximately  10  times  faster;  however,  the  sidelobe 
was  3  dB  worse. 

5.  Results  of  a  Lower-Contrast  Target  Case 

A.  Experimental  Results  of  a  Dense  Array 
To  study  the  effects  of  target  contrast  on  the  quality 
of  the  reconstructed  image  for  each  array  configura¬ 
tion,  we  conducted  a  set  of  experiments  with  a  lower- 
contrast  target  (|xa  =  0.10  cm-1)  embedded  in  the 
Intralipid.  The  center  of  the  target  was  located  at 
(*  =  0,  y  =  0,  2  =  2.5  cm).  The  centers  of  the 
imaging  voxels  in  z  were  0.5, 1.5,  2.5,  and  3.5  cm  for 


Fig  12  Experimental  images  of  ji0  at  target  layer  3  recon¬ 
structed  from  a  total  of  24  X  5  =  120  source- detector  pairs,  (a) 
Reconstructed  |ia  at  target  layer  3  with  510  iterations.  The  mea¬ 
sured  was  0.049  cm'1,  which  was  49%  of  the  true  target  |i0, 

and  its  location  was  (*  =  0.5,  y  =  -1.0),  which  was  displaced  from 
the  true  target  location  by  1.11  cm.  Image  artifacts  were  ob¬ 
served,  and  the  peak  was  -3  dB  from  the  (b)  -6-dB 

contour  plot  of  (a),  (c)  Reconstructed  |i0  at  target  layer  3  with 
1500  iterations.  The  peak  artifact  was  5  dB  higher  than  the 
image  lobe,  (d)  Reconstructed  |i„  at  target  layer  3  (56  iterations) 
with  only  a  prion  target  depth  information  provided  by  ultra¬ 
sound.  The  measured  |la(max)  was  0.074  cm'1,  and  its  Nation 
was  (x  =  0.5,  y  =  -0.5).  Image  artifacts  were  observed,  and  the 
peak  was  -5  dB  from  the  |ia(raaxv 


layers  1,  2,  3,  and  4,  respectively.  Figure  11  shows 
the  images  of  the  lower-contrast  target  obtained  at 
target  layer  3  [Fig.  11(a)]  and  deeper  nontarget  layer 
4  [Fig.  11(b)].  The  images  were  reconstructed  from 
the  measurements  made  with  the  dense  array.  At 
target  layer  3,  the  measured  |ia(ma3c) was  0063  c™,  ’ 
which  was  approximately  63%  of  the  target  p.a.  Ine 
measured  spatial  location  of  |ia(max) was  - 

-0.5),  which  agreed  reasonably  well  with  the  true 
target  location.  The  measured  -6-dB  beam  width 
was  1.83  cm,  which  was  144%  times  broader  than  the 
beam  width  of  the  high-contrast  case.  Edge  arti¬ 
facts  were  observed,  and  the  peak  level  was 
from  the  |ia(max).  At  nontarget  layer  4,  the  measured 

P-a(max)  waS  °-0871  CH1  >  whkh  W&S  eVen 

that  measured  at  the  target  layer.  Because  the  SNR 
of  the  data  was  lower  than  the  high-contrast  case,  the 
target  depth  estimate  was  poorer.  A  target  o 
,1  ,  ,  =  0.0543  cm"1  located  at  (x  =  0,  y  -  -O.h) 

wasTobserved  at  nontarget  layer  2,  and  its  mass  was 
much  smaller  than  that  obtained  at  the  target  layer 

and  nontarget  layer  4.  ,  ,  ,  , 

Figure  11(c)  is  the  reconstructed  pa  at  the  target 
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layer  from  only  the  dense  array  measurements.  A 
total  of  56  iterations  were  used  to  obtain  jiQ(max)  = 
0.107  cm-1,  and  the  reconstruction  was  approxi¬ 
mately  ten  times  faster  than  that  without  the  depth 
information.  The  spatial  resolution  was  8%  better 
than  that  obtained  from  Fig.  11(a),  and  the  -6-dB 
beam  width  was  1.67  cm.  The  contrast  resolution 
was  2  dB  worse. 

B.  Experimental  Results  of  a  Sparse  Array 

The  imaging  quality  of  sparse  arrays  decreased. 
Figure  12(a)  is  an  image  of  the  same  target  (pa  =  0.10 
cm  J)  reconstructed  from  measurements  made  with 
the  24  X  5  sparse  array.  The  measured  pa(max)  was 
0.049  cm  \  which  was  49%  of  the  true  value.  The 
-6-dB  contour  plot  is  shown  in  Fig.  12(b).  The  mea¬ 
sured  spatial  location  of  jia(max)  was  (x  =  0.5,  y  = 
-1.0),  which  was  displaced  from  the  true  target  lo¬ 
cation  by  1.11  cm  in  radius.  The  measured  -6-dB 
beam  width  was  2.98  cm,  which  was  163%  broader 
than  that  measured  from  the  dense  array.  Side- 
lobes  were  abundant,  and  the  peak  value  was  -3  dB 
below  the  peak  of  the  image  lobe.  These  sidelobes 
would  produce  false  targets  in  the  image  if  no  a  priori 
information  about  the  target  locations  were  given. 

In  this  case,  continuous  iteration  did  not  increase 
the  target  strength  but  increased  the  sidelobe 
strength.  Figure  12(a)  was  obtained  at  510  itera¬ 
tions,  whereas  Fig.  12(c)  was  obtained  at  1500  itera¬ 
tions.  After  approximately  three  times  more 
iterations,  the  peak  of  the  artifact  was  5  dB  higher 
than  the  peak  of  the  image  lobe.  The  target  depth 
estimated  from  3-D  images  was  worse  at  1500  itera¬ 
tions  than  that  at  the  510  iterations.  The  measured 
target  strengths  at  nontarget  layer  2  were  0.064  and 
0.1661  cm  1  at  510  and  1500  iterations,  respectively, 
and  the  strengths  at  nontarget  layer  4  were  0.0597 
and  0.1087  cm  \  respectively.  In  addition,  the 
background  noise  fluctuation  or  standard  deviation 
measured  at  nontarget  layer  1  was  increased  with  the 
iterations.  The  mean  and  the  standard  deviation  at 
510  iterations  were  0.024  and  0.0041  cm-1,  whereas 
these  values  at  1500  iterations  were  0.0267  and 
0.0104  cm  \ 

As  shown  in  Fig.  12(d),  the  reconstructed  image 
improved  a  lot  when  the  target  depth  was  given. 
The  maximum  strength  was  0.074  cm-1,  and  its  lo¬ 
cation  was  (0.5,  -0.5).  The  -6-dB  beam  width  was 
2.63  cm,  which  was  12%  better  than  that  obtained 
from  Fig.  12(a).  The  sidelobe  was  -6  dB  from  the 
peak,  which  was  improved  by  3  dB  compared  with 
Fig.  12(a). 

.  C.  Experimental  Results  of  Reconstructed  (la(max)  versus 
Total  Number  of  Source-Detector  Pairs 

With  the  iteration  number  normalized  to  the  dense 
array  case,  we  measured  target  strengths  at  target 
layer  3  and  nontarget  layer  4  for  all  sparse  array 
configurations.  Figure  13  shows  the  experimental 
data  points  (circles)  of  measured  |ia(max)  versus  the 
total  number  of  source- detector  pairs  obtained  at 
target  layer  3  [Fig.  13(a)]  and  nontarget  layer  4  [Fig. 
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Fig.  13.  Low-contrast  target  case,  (a)  Reconstructed  fiQ(maX)  ver¬ 
sus  total  source- detector  pairs  with  the  target  depth  available 
(stars)  and  the  curve-fitting  results  (upper  curve).  Reconstructed 
Ma(max)  versus  total  source— detector  pairs  measured  at  target  layer 
(circles)  and  the  curve-fitting  results  (lower  curve)  and  (b)  at 
deeper  nontarget  layer  4. 


13(b)],  The  two  curves  were  the  fitting  results  of 
experimental  points  when  we  used  second-order  poly¬ 
nomials.  In  both  layers,  the  |ia(max)  values  were  de¬ 
creased  when  the  total  number  of  source- detector 
pairs  was  reduced.  The  reconstructed  target 
strengths  were  reduced  to  less  than  60%  when  the 
total  pairs  were  less  than  156  and  140  for  target  layer 
3  and  nontarget  layer  4,  respectively.  The 
ultrasound-assisted  reconstruction  results  are  shown 
in  Fig.  13(a)  (stars),  and  the  average  reconstructed 
Ao(max)  for  all  array  configurations  was  0.085  cm"1. 
Compared  with  the  average  of  0.055  cm"1  obtained 
from  optical  imaging  only  at  the  target  layer,  a  30% 
improvement  was  achieved. 

D.  Experimental  Results  of  Imaging  Parameters  versus 
Total  Number  of  Source-Detector  Pairs 

The  measured  imaging  parameters  obtained  from  dif¬ 
ferent  array  configurations  are  listed  in  Table  2. 
Similar  to  Table  1,  listed  are  the  measured  parame¬ 
ters  at  target  layer  3  and  nontarget  layer  4.  The 
increase  in  beam  width  at  the  target  layer  was  neg¬ 
ligible  for  the  arrays  with  more  than  156  source- 
detector  pairs  and  was  more  than  50%  for  the  sparse 
an'ays  with  total  pairs  less  than  this  number.  The 
sidelobes  were  progressively  worse  when  the  total 


nairs  were  reduced.  At  the  target  layer,  the  distance 
between  the  |i0(maX)  location  and  the  true  target  loca¬ 
tion  was  more  than  1  cm  for  the  arrays  with  total 
pairs  less  than  140.  Table  2  then  lists  the  measured 
peak  image  lobe  at  nontarget  layer  2  and  its  location. 
The  target  appeared  at  nontarget  layer  2  in  all  array 
configurations  because  of  the  lower  SNR  of  the  data. 
However,  the  target  mass  observed  at  this  layer  was 
much  smaller  than  that  at  layers  3  and  4  in  all  cases. 
The  last  row  in  Table  2  shows  the  measured  imaging 
parameters  when  the  target  depth  was  available  to 
optical  reconstruction.  Compared  with  parameters 
obtained  from  optical  imaging  only,  the  -6-dB  beam 
width  was  improved  by  41%  on  average,  and  the  re¬ 
construction  speed  was  approximately  ten  times  fast¬ 
er;  however,  the  sidelobe  was  1  dB  worse. 

6.  Discussion 

In  addition  to  the  total  number  of  source- detector 
pairs,  the  measured  imaging  parameters  are  also  re¬ 
lated  to  other  system  parameters,  for  example,  mod¬ 
ulation  frequency  and  system  noise.  The  140-MHz 
modulation  frequency  chosen  in  this  study  is  a  typical 
frequency  used  by  many  research  groups.  Our  sys¬ 
tem  noise  level,  including  both  coherent  and  incoher¬ 
ent,  is  less  than  10  mV  peak  to  peak,  which  is 
sufficiently  low.  Therefore  the  results  we  obtained 
are  pertinent  to  3-D  imaging  using  similar  system 
parameters  and  reflection  geometry. 

In  this  study  the  target  absorption  coefficient  was 
reconstructed  from  measurements.  Similar  studies 
can  be  done  for  the  scattering  coefficient  as  well.  To 
reconstruct  the  scattering  coefficient,  we  can  use  the 
scattering  weight  matrix  derived  by  O’Leary11  to  re¬ 
late  the  medium  scattering  variations  with  the  mea¬ 
surements.  Simultaneous  reconstruction  of  both 
absorption  and  scattering  coefficients  in  reflection  ge¬ 
ometry  is  also  possible,  provided  that  the  absorption 
and  scattering  weight  matrices  are  regulated  care¬ 
fully.  Because  the  eigenvalues  of  the  two  matrices 
are  significantly  different,  good  regulation  schemes 
are  needed  to  balance  the  reconstructed  absorption 
and  scattering  coefficients  at  each  iteration.  This 
subject  is  one  of  our  topics  for  further  study. 

In  this  study  the  ultrasound-assisted  optical  recon¬ 
struction  was  demonstrated  at  a  particular  target 
layer.  Similar  studies  can  be  done  with  multiple 
targets  located  at  different  layers.  In  the  multiple 
target  case,  we  can  attribute  measured  perturbations 
to  more  layers  instead  of  a  single  layer.  However, 
the  improvements  in  reconstructed  optical  properties 
and  reconstruction  speed  may  be  less  than  that  of  the 
single-layer  case. 

Ultrasound  has  good  imaging  capability,  and  it  can 
detect  small  lesions  of  a  few  millimeters  in  size. 
However,  its  specificity  in  cancer  detection  is  not  high 
as  a  result  of  overlapping  characteristics  of  benign 
and  malignant  lesions.  NIR  imaging  has  high  spec¬ 
ificity  in  cancer  detection;  however,  it  suffers  low  res¬ 
olution  and  lesion  location  uncertainty  because  of  the 
diffused  nature  of  the  NIR  light.  The  hybrid  imag¬ 
ing  that  combines  ultrasound  imaging  capability  and 


NIR  contrast  has  a  great  potential  to  overcome  defi¬ 
ciencies  of  either  method.  As  we  reported  in  the 
paper,  the  target  depth  information  can  significantly 
improve  the  accuracy  of  the  reconstructed  optical  ab¬ 
sorption  coefficient  and  reconstruction  speed.  In  ad¬ 
dition  to  use  of  a  priori  target  depth  information,  the 
target  spatial  distribution  provided  by  ultrasound 
can  be  used  in  optical  reconstruction  as  well.18 
With  the  localized  spatial  and  temporal  target  in¬ 
formation,  the  accuracy  of  the  reconstructed  optical 
properties  and  the  reconstruction  speed  can  be  im¬ 
proved  further.  To  demonstrate  this,  we  need  an 
ultrasound  imaging  transducer  located  at  the  mid¬ 
dle  portion  of  the  probe.  We  are  currently  pursu¬ 
ing  this  study.23-24 

In  this  study  the  targets  of  different  contrasts  were 
located  at  the  center  position.  We  have  also  done 
studies  with  targets  of  different  contrasts  located  at 
off-center  positions.  For  an  off-center  target  case, 
the  effective  number  of  source- detector  pairs  is  less 
than  that  of  the  on-center  target  case  because  mea¬ 
surements  from  certain  source— detector  pairs  do  not 
contain  much  information  about  the  target.  For  ex¬ 
ample,  if  a  target  is  placed  at  (x  =  2,  y  =  2,  z  =  3.0 
cm),  the  measurements  of  source- detector  pairs  at 
the  opposite  comer  of  the  probe  contribute  little  to  the 
image  reconstruction.  In  one  study,  targets  of  high 
and  low  contrast  were  located  at  (x  =  2,y  =  2,  z  =  3.0 
cm).  The  reconstructed  maximum  absorption  coef¬ 
ficient  at  the  target  layer  was  related  more  to  the 
total  neighbor  source- detector  pairs.  However,  the 
maximum  value  at  a  deeper  nontarget  layer  was  re¬ 
lated  more  to  the  total  source- detector  pairs  and  was 
decreased  with  the  reduction  of  total  pairs.  Because 
the  photons  originated  from  distant  sources  and  de¬ 
tected  by  distant  detectors  experience  longer  and 
more  diffused  scattering  paths,  they  are  likely  to  in¬ 
teract  with  the  off-center  target  and  contribute  to  the 
absorption  estimate  at  the  deeper  layer.  In  the 
same  study,  the  measured  sidelobe  levels  were  4.5 
and  2.0  dB  poorer  on  average  compared  with  the 
on-center  high-  and  low-target  cases,  respectively. 
The  beam  widths  were  comparable  to  those  measured 
from  on-center  cases. 

7.  Summary 

The  relationship  between  the  total  number  of  source- 
detector  pairs  and  the  imaging  parameters  of  a  re¬ 
constructed  absorption  coefficient  was  evaluated 
experimentally.  A  frequency-domain  system  of  a 
140-MHz  modulation  frequency  was  used  in  the  ex¬ 
periments.  Reconstruction  at  a  selected  target 
depth  with  a  priori  depth  information  provided  by 
ultrasound  was  demonstrated.  The  results  have 
shown  that  the  reconstructed  absorption  coefficient 
and  the  spatial  resolution  of  the  absorption  image 
were  decreased  when  the  total  number  of  source- 
detector  pairs  was  reduced.  More  than  160  source- 
detector  pairs  were  needed  to  reconstruct  the 
absorption  coefficient  within  60%  of  the  true  value 
and  spatial  resolution  comparable  to  that  obtained 
with  the  dense  array.  The  contrast  resolution  was 
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'  poorer  in  general  because  of  edge  artifacts  and  could 
be  worse  if  significant  larger  iteration  numbers  are 
used  for  reconstruction.  The  error  in  target  depth 
estimated  from  3-D  optical  images  was  approxi¬ 
mately  1  cm.  With  the  a  priori  target  depth  infor¬ 
mation  provided  by  ultrasound,  the  reconstruction 
can  be  done  at  a  selected  depth.  Because  the  un¬ 
knowns  were  reduced  significantly,  the  reconstruc¬ 
tion  speed  was  approximately  ten  times  faster  than 
that  without  depth  information.  In  addition,  the  ac¬ 
curacy  of  the  reconstructed  absorption  coefficient  was 
improved  by  15%  and  30%  on  average  for  high-  and 
low-contrast  cases,  respectively.  Furthermore,  the 
measured  -6-dB  beam  width  was  improved  by  24% 
and  41%  for  high-  and  low-contrast  cases,  respec¬ 
tively.  The  sidelobe  was  3  and  1  dB  poorer  for  high- 
and  low-contrast  cases  because  the  measurement 
noise  was  lumped  to  single-layer  reconstruction  in¬ 
stead  of  multiple  layers. 

In  conclusion,  ultrasound-assisted  3-D  optical  im¬ 
aging  has  shown  promising  results  to  overcome  the 
problems  associated  with  the  reconstruction  by  use  of 
diffusive  waves.  With  the  target  depth  information 
provided  by  ultrasound,  the  reconstructed  absorption 
coefficient  was  more  accurate  and  the  reconstruction 
speed  was  much  faster. 
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